# Measuring correlations from the collective spin fluctuations of a large
ensemble of lattice-trapped dipolar spin-3 atoms

Youssef Aziz Alaoui1,2, Bihui Zhu3, Sean Robert Muleady4,5, William
Dubosclard1,2, Tommaso Roscilde6, Ana Maria Rey4,5, Bruno Laburthe-Tolra2,1,
and Laurent Vernac1,2 1 Université Paris 13, Laboratoire de Physique des
Lasers, F-93430, Villetaneuse, France
2 CNRS, UMR 7538, LPL, F-93430, Villetaneuse, France
3 Homer L. Dodge Department of Physics and Astronomy, The University of
Oklahoma, Norman, Oklahoma 73019, USA &
Center for Quantum Research and Technology, The University of Oklahoma,
Norman, Oklahoma 73019, USA
4 JILA, NIST and Department of Physics, University of Colorado, Boulder, USA
5 Center for Theory of Quantum Matter, University of Colorado, Boulder, CO
80309, USA
6 Univ Lyon, Ens de Lyon, CNRS, Laboratoire de Physique, F-69342 Lyon, France

###### Abstract

We perform collective spin measurements to study the buildup of two-body
correlations between $\approx 10^{4}$ spin $s=3$ chromium atoms pinned in a 3D
optical lattice. The spins interact via long range and anisotropic dipolar
interactions. From the fluctuations of total magnetization, measured at the
standard quantum limit, we estimate the dynamical growth of the connected
pairwise correlations associated with magnetization. The quantum nature of the
correlations is assessed by comparisons with short and long time expansions,
and numerical simulations. Our work shows that measuring fluctuations of spin
populations provides new ways to characterize correlations in quantum many-
body systems, for $s>1/2$ spins.

Characterizing quantum correlations between different parts of a system is of
fundamental importance for the development of quantum technologies and the
study of complex quantum systems. Quantum correlations are at the heart of the
most peculiar effects predicted by quantum mechanics, such as entanglement,
EPR steering Horodecki  _et al._ (2009); Cavalcanti and Skrzypczyk (2016);
Pezzè  _et al._ (2018), or Bell nonlocality Brunner  _et al._ (2014), which
all give advantage for different quantum information or metrological tasks.
Quantum correlations can even arise for non-entangled states Ollivier and
Zurek (2001); Luo (2008), where they can still constitute an interesting
resource Bera  _et al._ (2017). Systems made of $s>1/2$ particles pinned in
optical lattices are particularly interesting for such applications, as their
Hilbert space, enlarged with respect to qubit ($s=1/2$) systems, offers new
possibilities for quantum information processing Wang  _et al._ (2020).

Quantum correlations should appear in generic quantum systems Ferraro  _et
al._ (2010), but proving their inherent quantum nature is an experimental
challenge, which requires the measurement of non-commuting operators. As full
state-tomography scales exponentially with the number of constituents Kaufman
_et al._ (2016) and thus becomes impossible in large ensembles, it is of
crucial importance to develop new protocols to infer correlations from partial
measurements such as bipartite or collective measurements. The latter have
been successful in demonstrating entanglement Pezzè  _et al._ (2018),
steeringFadel  _et al._ (2018); Kunkel  _et al._ (2018); Lange  _et al._
(2018), or nonlocality Schmied  _et al._ (2016), in experimental platforms
dealing with effective two-level systems, for which entanglement witnesses
have a simpler structure compared to particles with $s>1/2$ Sørensen and
Mølmer (2001); Vitagliano  _et al._ (2011, 2014). Extensions to $s=1$ systems
in spinor Bose-Einstein condensates have demonstrated number squeezing in pair
creation processes via spin-mixing collisions Pezzè  _et al._ (2018); Bookjans
_et al._ (2011); Gross  _et al._ (2011); Lücke  _et al._ (2011); Qu  _et al._
(2020), SU(1,1) interferometryLinnemann  _et al._ (2016), and entangled
fragmented phases Evrard  _et al._ (2021). These systems nevertheless operated
in the regime where the single-mode approximation is valid Law  _et al._
(1998), which enormously simplifies the quantum dynamics.

Figure 1: Principle of the experiment. a) The cartoons zoom over a small
region of the Mott insulating distribution with doubly (core) and simply
occupied sites (shell). Spin 3 chromium atoms are excited at $t=0$ by a RF
pulse, with 5 cycles at the Larmor period set by the external magnetic field
B. The spin directions then makes an angle $\theta$ (set to $\pi/2$ in this
work) with respect to B, which triggers spin dynamics. Correlations develop
between spins, while doubly occupied sites get empty. b) Stern Gerlach
separation provides measurement of the fractional spin populations
$p_{ms}(t)$, through fluorescence imaging, at a given dynamic time $t$ .
Repeating the experiment allows us to compute the variance of the
magnetization, and hence the correlator $C_{z}$ of Eq.(1).

In this work we measure for the first time two-body correlations in a
macroscopic array of spin-3 chromium atoms pinned in a 3D optical lattice and
coupled via long-range and anisotropic magnetic dipolar interactions. Prior
experiments measuring out of equilibrium spin dynamics in these arrays
demonstrated compatibility with the growth of quantum correlations Lepoutre
_et al._ (2019); Patscheider  _et al._ (2020) and their approach to quantum
thermalization Lepoutre  _et al._ (2019). Here, we make use of the large
atomic spin to obtain a direct measurement of two-body correlations.
Specifically, after triggering out-of equilibrium spin dynamics, we acquire
statistics on the $2s+1=7$ spin populations, and quantify the growth of inter-
atomic spin correlations by analyzing the statistical fluctuations of the
collective spin component along the external magnetic field, i.e. the
magnetization. The quantum nature of the correlations we measure is validated
by agreement with exact short time expansions, with a high temperature series
expansion applied to the asymptotic quantum thermalized state at long time,
and with simulations of the full quantum dynamics via advanced phase-space
numerical methods.

We consider a system of $N$ spin $s>1/2$ particles. We define
$\hat{s}_{z}^{i}$ as the z-component of the spin of the $i^{th}$ particle. The
correlator $C_{z}$ we aim to measure is:

$C_{z}=\sum_{i\neq
j}^{N}\left(\left\langle\hat{s}_{z}^{i}\hat{s}_{z}^{j}\right\rangle-\left\langle\hat{s}_{z}^{i}\right\rangle\left\langle\hat{s}_{z}^{j}\right\rangle\right)={\rm{Var}}(\hat{S}_{z})-\Sigma_{z}$
(1)

with
${\rm{Var}}(\hat{S}_{z})=\langle\hat{S}_{z}^{2}\rangle-\langle\hat{S}_{z}\rangle^{2}$
the variance of the collective spin component
$\hat{S}_{z}=\sum_{i=1}^{N}\hat{s}_{z}^{i}$, and $\Sigma_{z}$ the sum of
individual variances,
$\Sigma_{z}=\sum_{i=1}^{N}\Big{(}\left\langle(\hat{s}_{z}^{i})^{2}\right\rangle-\left\langle\hat{s}_{z}^{i}\right\rangle^{2}\Big{)}$.
$\Sigma_{z}$ accounts for intraparticle correlations, which are only non-
trivial for $s>1/2$, as $\left\langle({\hat{s}}_{z}^{i})^{2}\right\rangle=1/4$
if $s=1/2$. The interparticle correlations are accounted for by the two-body
correlator $C_{z}$.

In this work, we independently determine ${\rm{Var}}(\hat{S}_{z})$ and
$\Sigma_{z}$ from collective measurements in order to show that $C_{z}$
departs from 0. Measurement of ${\rm Var}(\hat{S}_{z})$ requires the
experiment to be repeated many times to acquire adequate statistics.
Measurement of $\Sigma_{z}$ is straightforward in the case of an homogeneous
system, comprising singly-occupied lattice sites, referred as singlons in the
following. Indeed for singlons,
$\left\langle(\hat{s}_{z}^{i})^{2}\right\rangle=\sum_{m_{s}}p_{m_{s}}^{(i)}m_{s}^{2}$,
with $p_{m_{s}}^{(i)}$ the probability that the site $i$, uniquely populated
by the $i$-th spin, is in the $m_{s}$ spin state ($-3\leq m_{s}\leq 3$ ,
$\sum_{m_{s}}p_{m_{s}}^{(i)}=1$) so that
$\sum_{i}\left\langle(\hat{s}_{z}^{i})^{2}\right\rangle=N\sum_{m_{s}}p_{m_{s}}m_{s}^{2}$
even without homogeneity. Homogeneity ensures that $p_{m_{s}}^{(i)}=p_{m_{s}}$
are independent of site $i$, so that
$\left\langle\hat{s}_{z}^{i}\right\rangle^{2}=(\sum_{m_{s}}p_{m_{s}}m_{s})^{2}$;
therefore

$\frac{\Sigma_{z}}{N}=\sum_{m_{s}}p_{m_{s}}m_{s}^{2}-\left(\sum_{m_{s}}p_{m_{s}}m_{s}\right)^{2}$
(2)

We show in Ref. Sup  that the experimental magnetic inhomogeneities lead to
negligible deviations from Eq.(2). Besides, inhomogeneities of the lattice
potential are below $2.5\%$, and therefore have a negligible effect on spin
dynamics in the Mott regime Fersterer  _et al._ (2019). Finally, we show in
Ref. Sup  that Eq.(2) still holds for doubly occupied sites (doublons) which
are populated at the beginning of dynamics (see below). This is due to the
fact that the spin of each pair of particles in a doubly-occupied site is
well-defined at all times. Therefore, $\Sigma_{z}$ is given by Eq.(2) for the
whole dynamics, and measurement of $p_{m_{s}}=N_{m_{s}}/N$, with $N_{m_{s}}$
the total number of atoms in spin state $m_{s}$, yields $\Sigma_{z}$.

Experimental setup. The starting point of our experiments is a spin-polarized
52Cr Bose-Einstein Condensate (BEC) produced in a crossed dipole trap, with
typically $15000$ atoms polarized in the minimal Zeeman energy state
$m_{s}=-3$. We load the 52Cr BEC in a 3D optical lattice deep into the Mott
insulator regime. The lattice implemented with five lasers at
$\lambda_{L}=532$ nm is described in Lepoutre  _et al._ (2019). The total
lattice depth is equal to 60 recoil energy at $\lambda_{L}$. We estimate the
tunneling time to be $\simeq 20$ ms. We obtain a core of doublons comprising
$\simeq 50\%$ of the atoms, surrounded by a shell of singlons.

As shown in Fig. 1(a), we trigger spin dynamics by rotating all spins with the
use of a Radio Frequency (RF) $\pi/2$ pulse. After the pulse all spins are
oriented orthogonal to the external magnetic field. The Larmor frequency
$f_{L}=g_{L}\mu_{B}B_{0}/\hbar$ (with $g_{L}\simeq 2$ the Landé factor,
$\mu_{B}$ the Bohr magneton, and $B_{0}=0.75$ Gauss the amplitude of the
magnetic field) is $f_{L}\simeq 2.1$ MHz. The RF frequency $f_{RF}$ is set at
resonance, and fluctuations of the detuning $(f_{L}-f_{RF})\simeq 1$ kHz are
small compared to the RF Rabi frequency $f_{R}$, thanks to the use of a 30
Watt RF amplifier. In practice the RF pulse has a duration of exactly 5 Larmor
periods, with $f_{R}=\frac{1}{5}\frac{f_{L}}{4}=105$ kHz; the
$\bar{\theta}=\pi/2$ pulse is set to have an identical initial phase at each
realization. Fluctuations of the rotation angle $\theta$ are estimated to have
a standard deviation of $\sigma_{\theta}\simeq 2.5\times 10^{-3}$ rad (see
below). After the initial state preparation with the RF pulse, spins interact
via magnetic dipolar interactions in the optical lattice for a duration $t$.
We then adiabatically ramp down the optical lattice, and proceed to
measurements.

Theoretical model Dipolar interactions between singlons during the dark time
evolution are described by the effective dipolar Hamiltonian $\hat{H}_{\rm
dd}$, which is a XXZ spin model Hamiltonian:

$\hat{H}_{\rm
dd}=\sum_{i>j}^{N}V_{ij}\left[\hat{s}_{z}^{i}\hat{s}_{z}^{j}-\frac{1}{2}\left(\hat{s}_{x}^{i}\hat{s}_{x}^{j}+\hat{s}_{y}^{i}\hat{s}_{y}^{j}\right)\right]$
(3)

with $V_{ij}=V_{dd}\left(\frac{1-3\cos^{2}\theta_{ij}}{r_{ij}^{3}}\right)$,
$V_{dd}=\frac{\mu_{0}(g_{L}\mu_{B})^{2}}{4\pi}$, and $\mu_{0}$ the magnetic
permeability of vacuum. The sum runs over all pairs of particles ($i$,$j$),
$r_{ij}$ is their corresponding distance, $\theta_{ij}$ the angle between
their inter-atomic axis and the external magnetic field,
${\bf{\hat{s}}}_{i}=\\{\hat{s}_{x}^{i},\hat{s}_{y}^{i},\hat{s}_{z}^{i}\\}$ are
spin-3 angular momentum operators for atom $i$. The shortest intersite
distance $r_{min}=268$ nm in our lattice Lepoutre  _et al._ (2019) corresponds
to a dipolar coupling $V_{dd}/r_{min}^{3}\simeq h\times 3$ Hz.

Figure 2: Evolution of the atom number $N(t)$ (Top), and of the fluctuations
measured in the experiment (Bottom): we show the standard deviations of the
normalized magnetization $M_{z}$ and of the technical noises featured in
Eq.(4): $\sigma_{{\rm{prep}}}$ (shaded area; preparation of the sample),
$\sigma_{{\rm{fit}}}$ (fitting uncertainties), and $\sigma_{{\rm{shotnoise}}}$
(fluorescence imaging). The quantum projection noise
$\sqrt{\frac{3}{2N(t)}}=\sigma_{{\rm{SQN}}}$ is shown for comparison. Lines
are guides to the eye. Error bars evaluated from statistics correspond to two
standard deviations.

Given the strong contact interactions that favor spin alignment Lepoutre  _et
al._ (2018a, b) and the fully polarized initial state, the same Hamiltonian
can be used to describe the dynamics of doubly occupied sites (doublons) just
by replacing ${\bf{\hat{s}}}_{i}$ by a spin-6 angular momentum operator at the
corresponding site Fersterer  _et al._ (2019). Nevertheless, as soon as the
spin excitation is performed, doublons start to leave the trap, due to dipolar
relaxation de Paz  _et al._ (2013), see Fig.2: for $0<t<10$ ms the spin system
comprises both singlons and doublons, but only singlons remain for $t>10$ ms
and losses become negligible. This is why we run simulations for singlons
only, which allows for quantitative comparison with the experiment except at
short time.

As shown in previous work Lepoutre  _et al._ (2019), we need to include the
one-body term $\hat{H}_{Q}=B_{\rm Q}\sum_{i}^{N}(\hat{s}_{z}^{i})^{2}$
accounting for light shifts created by the lattice lasers. Finally, spin
dynamics is also driven at some point by tunneling processes. However in the
Mott regime tunneling-assisted superexchange processes are happening at longer
time scales and remain irrelevant for the current measurements.

Correlation dynamics During the evolution under $\hat{H}_{\rm
dd}+\hat{H}_{Q}$, $\langle\hat{S}_{z}\rangle$ and ${\rm{Var}}(\hat{S}_{z})$
are constant, as these two operators commute with $\hat{S}_{z}$; on the
contrary, interactions between spins lead to evolution of spin populations,
hence of $\Sigma_{z}$ and $C_{z}$. In our case, as spins are orthogonal to the
magnetic field, $\langle\hat{S}_{z}\rangle=0$ and Eq.(2) reads
$\Sigma_{z}=N\sum_{m_{s}}p_{m_{s}}m_{s}^{2}$; besides
${\rm{Var}}(\hat{S}_{z})=\frac{3}{2}N$ as the initial state is a coherent spin
state. The short-time evolution is obtained by perturbation theory Lepoutre
_et al._ (2019), leading to $C_{z}\approx-\frac{45N}{8}t^{2}(3V_{\rm
eff}^{2}-4B_{Q})$, where $V_{\rm eff}^{2}=\sum_{i,j\neq
i}^{N}V_{ij}^{2}/(2N)$, $V_{\rm eff}\simeq h\times 4.3$ Hz. At longer times,
we can numerically simulate the dynamics via a semiclassical phase space
method known as the generalized discrete truncated Wigner approximation
(GDTWA) Zhu  _et al._ (2019), which was previously shown to capture
quantitatively the spin population dynamics of this system Lepoutre  _et al._
(2019).

We also provide a theoretical estimate of the expected correlation at long
time assuming the Eigenstate Thermalization Hypothesis D’Alessio  _et al._
(2016); Kaufman  _et al._ (2016). In this case, due to the build up of quantum
correlations, local observables at long time can be described by a thermal
density matrix with additional Lagrange multipliers that account for conserved
quantities. A high-temperature $T$ series expansion valid for our system
Lepoutre  _et al._ (2019); Bihui Zhu and Rey  leads to
$C_{z}(t\to\infty)=\left(-\frac{5}{2}+12\beta B_{\rm Q}\right)N$, with
$\beta\equiv\frac{1}{k_{B}T}=\frac{5B_{\rm Q}+9\bar{V}}{48V_{\rm
eff}^{2}+24B_{\rm Q}^{2}}$, and $\bar{V}=\frac{1}{N}\sum_{i>j}^{N}V_{ij}\simeq
h\times-0.6$ Hz.

Experimentally, the quantities of interest are the total number of atoms,
$N(t)$, and the fractional spin populations $p_{m_{s}}(t)$. While the
fluctuations in $N(t)$ from shot to shot (with a standard deviation of about
$10\%$) yield large extra fluctuations on the measured absolute spin
populations $N_{m_{s}}=Np_{m_{s}}$, this extra source of noise is cancelled
when dealing with fractional populations. For measuring the total atom number,
we use absorption imaging of the BEC. We checked that the loading in the
optical lattice does not lead to losses and therefore $N(t=0)$ is equal to the
atom number in the BEC. We estimate the accuracy of this measurement equal to
$10\%$.

To measure $p_{m_{s}}(t)$, we spatially separate the $7$ spin components
during a time of flight of $14$ ms, using a Stern Gerlach (SG) technique. We
use fluorescence imaging to count atoms: it brings equal detectivity of spin
components, and makes the use of Electron Multiplying (EM) on the CCD camera
advantageous for signal to noise ratio (see Sup  for details). Atoms are
excited by a saturating laser set at 425 nm (with a transition rate
$\Gamma=2\pi\times 5\times 10^{6}$ Hz) during typically $500$ $\mu$s. The
magnetic field $B_{0}$ is reduced to a small value ($g\mu_{B}B_{0}\ll
h\Gamma$) to ensure that the fluorescence rates of the 7 spin components are
almost equal. We use a “delta-kick” stage Ammann and Christensen (1997) at the
very beginning of the time of flight: it consists of a short 0.5 ms pulse of
an intense IR laser along the separation axis of the SG that applies a force
on the atoms and helps reducing velocity dispersion. We fine-tune the
frequency of the laser exciting the atoms and the amplitude of all three
components of the magnetic field during the fluorescence stage. The obtained
regular shape of clouds (see Fig.1(b)) favors efficient fitting.

Figure 3: (a) Symbols are experimental values for the two contributions to the
correlator $C_{z}$ (see eq.(1)) normalized to atom number. Full line is
results of simulations for $\Sigma_{z}$, while the dashed line shows the
expected value for ${\rm{Var}}(\hat{S}_{z})$ (see text). (b) Experimental
value of the correlator $C_{z}$ normalized to the atom number (symbols), with
comparison to simulations (full line), and short time expansion (dotted line).
The dashed line corresponds to the calculated value in the quantum thermalized
state. Error bars evaluated from statistics correspond to two standard
deviations.

By fitting of the atomic clouds with a Gaussian function we obtain the values
of the number of counts $C_{m_{s}}$ detected for every spin components $m_{s}$
which set the value of $p_{m_{s}}=C_{m_{s}}/\sum_{m_{s}}C_{m_{s}}$. During the
dynamics, $N(t)$ is deduced by multiplying $N(0)$ by the ratio of the total
number of counts at $t$ and at $t=0$.

As explained above, ${\rm{Var}}(\hat{S}_{z})(t)$ is expected to be equal to
$\frac{3}{2}N(t)$ for a dipolar system without losses. But we do not assume
that this equality holds, we measure ${\rm{Var}}(\hat{S}_{z})(t)$ by thorough
investigations of all different sources of noise. In practice, we measure the
variance of the normalized magnetization of the sample,
$M_{z}=\sum_{m_{s}}p_{m_{s}}m_{s}$, $-3\leq M_{z}\leq 3$, from 3 to 5 sets of
40 pictures. In absence of noise,
${\rm{Var}}(\hat{S}_{z})=N\times{\rm{Var}}(M_{z})$. But at $t=0$, we obtain
$N(0)\times{\rm{Var}}(M_{z})\simeq 2\times\frac{3}{2}$, which shows that noise
processes come into play in our measurement of $M_{z}$: a proper determination
of ${{\rm Var}(\hat{S}_{z})}$ requires to evaluate their contribution
independently.

The noises on $M_{z}$ originate from fluctuations in the preparation angle
$\theta$, in the detection process (due to the Poissonian nature of light),
and in the evaluation of counts on the camera (related to error in the fitting
procedure). We denote their respective contribution to the standard deviation
on $M_{z}$ as $\sigma_{{\rm{prep}}}$, $\sigma_{{\rm{shotnoise}}}$ and
$\sigma_{{\rm{fit}}}$. These different noises are statistically independent,
so that

${\rm{Var}}(M_{z}(t))=\frac{{\rm{Var}}(\hat{S}_{z}(t))}{N^{2}(t)}+\sigma_{{\rm{shotnoise}}}^{2}(t)+\sigma_{{\rm{fit}}}^{2}(t)+\sigma_{{\rm{prep}}}^{2}$
(4)

from which we derive ${\rm{Var}}(\hat{S}_{z})(t)$ at any time $t$.

From a first principle calculation we can determine
$\sigma_{{\rm{shotnoise}}}$ from the average counts $C_{m_{s}}$ and camera
parameters, see details in Sup . Similarly, $\sigma_{{\rm{fit}}}$ is well
evaluated from data analysis. We use measurements at $t=0$ to evaluate the
last contribution, $\sigma_{{\rm{prep}}}$. Indeed, the initial sample
corresponds to an uncorrelated spin coherent state, for which
${\rm{Var}}(\hat{S}_{z})=\frac{3}{2}N(0)$ is guaranteed. The conservation of
magnetization during the whole spin dynamics ensures that
$\sigma_{{\rm{prep}}}$ is constant, as discussed in Sup ; we stress that the
contribution of the preparation noise becomes negligible at long time, see
Fig. 2.

The noise contributions as dynamics proceeds are shown in Fig. 2, and compared
to the one of atomic projection noise,
$\sigma_{\rm{SQN}}=\sqrt{\frac{3}{2N(t)}}$. We obtain
$\sigma_{{\rm{prep}}}=0.008\simeq 0.7\sigma_{\rm{SQN}}$. As
$\sigma_{{\rm{prep}}}$ scales like $s\times g_{L}$ and is independent of
$N_{0}$, while $\sigma_{\rm{SQN}}$ scales like $\sqrt{\frac{s}{N(0)}}$ we
stress the difficulty to get such a low value with large $N(0)$, large spin
$s=3$ and large Landé factor $g_{L}=2$. Figure 2 shows that
$\sigma_{{\rm{shotnoise}}}$ scales as $1/\sqrt{N(t)}$ as expected, and that
$\sigma_{{\rm{fit}}}$ has about the same scaling.

We show our measurements of ${\rm{Var}}(\hat{S}_{z})(t)/N(t)$ in Fig. 3(a).
The scatter of the data points around $\frac{3}{2}$ is comparable to the
average error bars for the different points. This indicates that deviations
compared to $\frac{3}{2}$ are not statistically significant throughout the
curve. This arises because the number of shots taken to estimate
${\rm{Var}}(\hat{S}_{z})$ at each time ($\simeq 150$) is not large enough for
the noise associated with finite data sampling to be negligible - a well-known
difficulty when estimating correlations from noise analysis.

As we measure a substantial growth for $\Sigma_{z}(t)/N(t)$, we can assert
that the correlator $C_{z}(t)$ significantly differs from zero for $t>20$ ms,
as directly shown in Fig.3(b). Fig. 3(a) shows a good quantitative agreement
between the measured $\Sigma_{z}(t)/N(t)$ and predictions from our GDTWA
simulations assuming only singlons, while Fig. 3(b) shows qualitative
agreement for the measured $C_{z}(t)$ with our short time expansion. The value
of the quadratic term $B_{Q}$ in simulations, $B_{Q}=h\times-5.1$ Hz, is
inferred from population analysis during the whole dynamics Sup ; it leads to
$\frac{C_{z}}{N}(t\to\infty)\simeq-1.3$, in good agreement with the data.

Our measurements thus quantify the amount of two-body correlations in the
expected highly correlated state reached at long time. Assuming translational
invariance and isotropic correlations decaying exponentially with a
correlation length $\xi$, the measured $C_{z}$ and $\Sigma_{z}$ at long time
can be related to the onset of correlations with $\xi\approx 0.3$ (in units of
the lattice spacing) Sup . This estimate represents a lower bound to the
actual correlation length (assuming concentration of correlations at short
distance); its rather small value is nonetheless compatible with the scenario
of thermalization at high-temperature.

We now discuss the influence of losses. As dipolar spin exchange dynamics
proceeds, doublons can get correlated with surrounding singlons, resulting in
a modification of singlons fluctuations. Therefore, quantum fluctuations of
the sample, and consequently its quantum correlations, may differ from the
singlon-only case. Taking losses into account rigorously is difficult and
would require new theoretical models to be developed, which is beyond the
scope of this paper. We discuss simple arguments in Sup  to estimate the
contribution of losses on ${\rm{Var}}(\hat{S}_{z})(t)$, and predict small
corrections, at the 10 percent level. An improved experimental resolution
would be necessary to show deviation from a fully unitary system.

In conclusion, we have measured the growth of correlations in a large ensemble
of interacting spins by analyzing the fluctuations of the collective
magnetization. This achievement illustrates the new possibilities offered by
$s>1/2$ species, and represents an important step towards understanding the
complex quantum many-body dynamics in state-of-the-art simulators of quantum
magnetism.

Acknowledgements: We acknowledge careful review of this manuscript and useful
comments from Thomas Bilitewski and Lindsay Sonderhouse. The Villetaneuse
group acknowledges financial support from CNRS, Conseil Régional d’Ile-de-
France under Sirteq Agency, Agence Nationale de la Recherche (project
ANR-18-CE47-0004), and QuantERA ERA-NET (MAQS project). A.M.R is supported by
the AFOSR grant FA9550-18-1-0319, AFOSR MURI, by the DARPA DRINQs grant, the
ARO single investigator award W911NF-19-1-0210, the NSF PHY1820885, NSF JILA-
PFC PHY-1734006 grants, and by NIST.

##  I Supplemental material

###  I.1 Simulations supporting the small contribution of inhomogeneities

Here we quantify the corrections induced by the position dependent magnetic
field by explicitly accounting for it in GDTWA simulations. We use numerical
GDTWA simulations to justify the validity of the homogeneous approximation
used in the main text. Explicitly we add a term of the form $\sum_{\bf
i}B_{\bf i}{\hat{s}}_{z}^{\bf i}$ to our simulations with $B_{\bf
i}=-4.4i_{x}-5.1i_{y}+8.7i_{z}$Hz, and ${\bf i}=(i_{x},i_{y},i_{z})$ are
integers labeling the 3D lattice coordinates. The chosen homogeneity is
consistent with the expected value from experiment. As a result of the
spatially varying field, the spins precess at different rate and
$\langle\hat{s}_{z}^{i}\rangle\neq\langle\hat{S}^{z}\rangle/N$.

In Fig. 4, we show the average correlators $C_{z}$ and $\Sigma_{z}$ generated
by the inhomogeneity during the quantum dynamics. Whilethere is a nonzero net
spin projection along $z$, which causes the difference between results
obtained assuming homogeneous $\langle\hat{s}_{z}^{i}\rangle$ (solid lines)
and those with actual $\langle\hat{s}_{z}^{i}\rangle$ (circles), it is much
smaller than the value of the relevant correlators.

###  I.2 Calculation of $\Sigma_{z}$ parameter for doubly occupied sites

Here we focus on the case of doubly occupied sites and explain ways to
determine the relevant correlations considered in this work in terms of
observables that can be properly measured in the experiment. Particularly, we
want to find out how to connect the intra-spin correla tions $\Sigma_{z}$ to
spin populations.

For an ensemble of $N_{d}$ doublons, one atom can become correlated with the
other in the same site in a different way than with atoms in different sites,
and the correlations can be split accordingly:

$\displaystyle\langle(\hat{S}_{z})^{2}\rangle$
$\displaystyle=\sum_{i}\langle(\hat{s}^{i}_{z})^{2}\rangle+\sum_{i,l\neq
i}\langle\hat{s}^{i}_{z}\hat{s}^{l}_{z}\rangle$
$\displaystyle=\sum_{i}\langle(\hat{s}^{i1}_{z}+\hat{s}^{i2}_{z})^{2}\rangle+\sum_{i,l\neq
i}\langle\hat{s}^{i}_{z}\hat{s}^{l}_{z}\rangle={\Sigma}_{aa}+{\Sigma}_{ab}+\mathcal{C}_{z},$
(5)

where $i,l$ label different sites,
$\hat{s}^{i}_{z}=\hat{s}^{i1}_{z}+\hat{s}^{i2}_{z}$ accounting for two atoms
in site $i$,
$\Sigma_{aa}=\sum_{i}\langle(\hat{s}^{i1}_{z})^{2}+(\hat{s}^{i2}_{z})^{2}\rangle$,
and $\Sigma_{ab}=2\sum_{i}\langle\hat{s}^{i1}_{z}\hat{s}^{i2}_{z}\rangle$.

Figure 4: Effect of inhomogeneity. In the presence of magnetic field
gradients, comparing the numerical results of $C_{z}/N$ (orange) and
$\Sigma_{z}/N$ (blue) assuming homogeneous
$\langle\hat{s}_{z}^{i}\rangle=\sum_{m_{s}}m_{s}p_{m_{s}}$ (solid lines) and
without such assumptions (circles) shows that for our experimental parameters
the effect of magnetic field gradients on the local magnetization is
sufficiently small to be ignored in the analysis of the collective spin
variance.

Given that the doublons start fully polarized $S=6$ at the initial time, then
the main assumption is that there is a large energy gap opened by the contact
interactions that prevents demagnetization of the local spin (see Ref.
Fersterer  _et al._ (2019)). Under this assumption we can obtain the dynamics
of populations on the thirteen different components, $N_{M}$ with
$M\to\\{-6,-5,\dots,5,6\\}$, of an array of $S=6$ spins, and then relate the
measured spin populations $p_{m_{s}}$ with $m_{s}\to\\{-3,-2,\dots,2,3\\}$ to
$N_{M}$:

$\displaystyle 2N_{d}p_{m_{s}}$ $\displaystyle=$ $\displaystyle
2\sum_{M=-6}^{6}[C^{6M}_{3m3M-m}]^{2}N_{M},$ (6)

Here, $C^{6M}_{3m3M-m}$ is a Clebsch-Gordon coefficient. To avoid confusion,
in the following, we will use $s$ and $m$ symbols for spin-3 operators, and
$S$, $M$ for spin-6 operators.

Since there are two atoms in the same site, there are correlations of the form
$\Sigma_{aa}^{i}\equiv\langle(\hat{s}^{i1}_{z})^{2}+(\hat{s}^{i2}_{z})^{2}\rangle$
for each site $i$. At a given lattice site we can write the doublon’s
wavefunction as $\left|\psi\right>=\sum_{M}c_{M}\left|M\right>$. Using it, we
can obtain the following expression for the intra-spin correlations at each
site:

$\displaystyle\Sigma_{aa}^{i}$ $\displaystyle=$
$\displaystyle\langle(\hat{s}^{i1}_{z})^{2}+(\hat{s}^{i2}_{z})^{2}\rangle$ (7)
$\displaystyle=$
$\displaystyle\sum_{M,M^{\prime}}c_{M}^{*}c_{M^{\prime}}\langle
SM|(\hat{s}^{i1}_{z})^{2}+(\hat{s}^{i2}_{z})^{2}|SM^{\prime}\rangle$
$\displaystyle=$
$\displaystyle\sum_{M,M^{\prime}}c_{M}^{*}c_{M^{\prime}}\sum_{m1,m2,m1^{\prime},m2^{\prime}}C^{SM}_{jm1jm2}\langle
m_{1}m_{2}|(\hat{s}^{i1}_{z})^{2}$ $\displaystyle+$
$\displaystyle(\hat{s}^{i2}_{z})^{2}|m_{1}^{\prime}m_{2}^{\prime}\rangle
C^{SM^{\prime}}_{jm1^{\prime}jm2^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{M,M^{\prime}}c_{M}^{*}c_{M^{\prime}}\sum_{m1,m2}C^{SM}_{jm1jm2}(m_{1}^{2}+m_{2}^{2})C^{SM^{\prime}}_{jm1jm2}$
$\displaystyle=$
$\displaystyle\sum_{m1,m2}\sum_{M}(m_{1}^{2}+m_{2}^{2})|c_{M}|^{2}[C^{SM}_{jm1jm2}]^{2}$
$\displaystyle=$
$\displaystyle\sum_{m1}m_{1}^{2}\sum_{m2}\sum_{M}N^{i}_{M}[C^{SM}_{jm1jm2}]^{2}$
$\displaystyle+$
$\displaystyle\sum_{m2}m_{2}^{2}\sum_{m1}\sum_{M}N^{i}_{M}[C^{SM}_{jm1jm2}]^{2}$
$\displaystyle=$ $\displaystyle
2\sum_{m1}m_{1}^{2}\sum_{m2}\sum_{M}N^{i}_{M}[C^{SM}_{jm1jm2}]^{2}$
$\displaystyle=$ $\displaystyle\sum_{m_{s}}{m_{s}}^{2}p^{i}_{m_{s}},$

and $\sum_{i}\Sigma_{aa}^{i}/N_{d}=\sum_{m_{s}}m_{s}^{2}p_{m_{s}}$. Therefore
even for the $S=6$ manifold, such a relation remains valid, as in the $S=3$
case. For a homogeneous system we thus again obtain
$\Sigma_{z}=N\sum_{m_{s}}p_{m_{s}}m_{s}^{2}-\langle\hat{S}_{z}\rangle^{2}/N$
even when there are doublons present. With this equation we can obtain the
spin-spin correlations $C_{z}=\mathcal{C}_{z}+\Sigma_{ab}$ from experimental
population measurements even in the presence of doublons. While we don’t
separately consider $\Sigma_{ab}$ and $\mathcal{C}_{z}$ in the main text, in
Sec. I.9, we provide their dynamics from simulations and show that they can be
readily obtained from experimental data with the knowledge of $\Sigma_{aa}$.

###  I.3 Details on Imaging System

We collect fluorescence light with a 2 inches diameter, 20 cm focal length
achromat. The collection efficiency of fluorescence light is $\simeq 4\times
10^{-3}$. The imaging system after the lens collecting fluorescence, made of
three other achromat lenses ensuring a magnification of 1.4, allows to match
the size of the full image (the 7 atomic clouds) with the full size of the CCD
chip of our EMCCD camera. The camera is cooled at $-90^{\circ}C$, which makes
black-body radiation negligible. The quantum efficiency of detection is 0.82
at 425 nm. We use an EM average gain of 24 (which we measure, see section
below), and a binning equal to 2 (the counts of 4 adjacent pixels are added
together). Optical shielding all over the imaging path leads to 0.7 photon of
stray light per pixel in average.

###  I.4 Derivation of the shot noise contribution to the fluctuations of
magnetization

To account for all the possible sources of fluctuations, one needs to consider
the physical process that is at play to estimate the number of atoms, i.e.
fluorescence imaging. Atoms emit light, which is an inherently random process
characterized by a Poisson distribution. These photons are collected by a
camera with quantum efficiency $\eta$, and then the signal is first amplified
by a gain $G_{1}$ referred usually as pre-gain (this first gain is in fact
typically a compression $G_{1}<1$ to optimize the camera well depth to the
dynamical range of the digital to analog converter), and then amplified again
using an electron amplifier of gain $G_{2}$, and finally collected.

Let us define $N_{\nu}$ as the photon number that impinges on the camera. This
is a fluctuating variable characterized by by a variance $\Delta
N_{\nu}=<N_{\nu}>$. We define $N_{e}$ the electronic signal created by the
photons arriving at one pixel. One photon creates one electron with a quantum
efficiency $\eta$. It is important to remember that this process is stochastic
and should be treated as such. There are in fact two independent fluctuating
processes: the Poisson fluctuations of light, and the excitation of one
electron by one photon in a camera pixel. The variance associated with both
these processes need to be calculated, and added in quadrature (since the
fluctuations of these process are independent).

For exactly $N_{\nu}$ photons impinging on a camera pixel, the non-unity
quantum efficiency $\eta$ of the detector results in a variance

$\Delta N_{e}^{(1)}=(\eta-\eta^{2})N_{\nu}.$ (8)

In addition to this variance, we need to consider the variance of the photon
number, when neglecting the stochastic nature of the detector:

$\Delta N_{e}^{(2)}=\eta^{2}\Delta N_{\nu}$ (9)

Therefore the total variance in the electron signal is:

$\Delta N_{e}=\Delta N_{e}^{(1)}+\Delta N_{e}^{(2)}=\eta<N_{\nu}>=\eta\Delta
N_{\nu}$ (10)

We thus simply deduce that $\Delta N_{e}/<N_{e}>=\Delta N_{\nu}/<N_{\nu}>=1$,
independent of quantum efficiency. Note that $\sqrt{\Delta
N_{e}}/<N_{e}>=\frac{1}{\eta}\sqrt{\Delta N_{\nu}}/N_{\nu}$, which indeed
shows that the signal to noise may be degraded when $\eta<1$.

We consider now the effect of a deterministic gain $G$. We simply have
$S=G<N_{e}>$ and $\sigma_{S}=G\sigma_{N_{e}}$, with
$\sigma_{N_{e}}=\sqrt{\Delta N_{e}}$. Therefore
$\sigma_{S}/S=\sigma_{N_{e}}/<N_{e}>$. The signal to noise is not modified by
the gain. Keeping in mind that $\Delta N_{e}/<N_{e}>=1$, we finally find:

$\sigma_{S}=\sqrt{G}\sqrt{S}$ (11)

As a consequence, when obtaining a signal $S$ on the intensified camera, there
is an associated fundamental quantum noise with standard deviation
$\sqrt{GS}$. This treatment applies to the pre-gain of the camera, assumed
deterministic, which implies that 1 photon gives $G_{1}$ electrons.

For a camera with an electron multiplying process, one has to consider the
stochastic nature of the gain process. For a large number of amplification
stages, and a large value of the EM gain $G_{2}$, this leads to an extra
factor $\sqrt{2}$ in the standard deviation, see formula (8) of reference
Robbins and Hadwen (2003):

$\sigma_{S}=\sqrt{2G}\sqrt{S}$ (12)

Taking into account the two gains $G_{1}$ and $G_{2}$, we obtain finally:

$\sigma_{S}=\sqrt{2G_{1}G_{2}}\sqrt{S}$ (13)

This shows that the whole amplification line needs to be well known in order
to predict accurately the contribution of the shot-noise. For our camera
settings, $G_{1}=\frac{1}{4.01}$ is expected from manufacturer data. We
optimized the EM gain $G_{2}$ for our experiment, with a corresponding
expected value $G_{2}=30$. As EM gains are known to be sensitive to aging, we
made the following measurements to infer a reliable value of the overall
experimental gain $G_{exp}=G_{1}G_{2}$.

First, we illuminate the camera with a laser and compare the average number of
counts for two camera settings: 1=conventional, i.e. no EM gain, with an
expected value of $G_{1}=3$; 2=with EM gain. The comparison between the
average number of counts in the two experiments leads to a value of the
average EM gain $G_{2}=24\pm 1$. This simple procedure does not validate the
value of $G_{1}$ in the EM mode.

The second method gives a value of the effective gain for each single
experimental shot. It relies on the analysis of dimly illuminated regions in
each of these shots. From these so called dark regions which extend over
thousands of pixels ($\geq 5000$), one can extract the probability
distribution of the number of counts per pixel $\mathcal{C}(x)$. Fitting this
empirical distribution with the appropriate model -see fig.5 \- and assuming
$G_{1}=\frac{1}{4}$ (as per manufacturer’s specifications), one can conclude
to the value of the $G_{2}$ gain. The model is summarized within the following
Cauchy product

$\mathcal{C}(x)=\sum_{i=-i_{max}}^{i_{max}}\mathcal{N}_{0,\sigma}(i)\times\sum_{j=-2}^{1}\mathcal{F}(4(x-i)+j)$
(14)

Where:

  * •

$\mathcal{F}(x)=p_{0}\delta_{x}+\sum_{n_{p}}p_{n_{p}}\frac{x^{n_{p}-1}e^{-\frac{x}{g}}}{g^{n_{p}}(n_{p}-1)!}$
describes the amplification process of Poisson distributed input charges
($p_{n_{p}}=\lambda^{n_{p}}\frac{e^{-\lambda}}{n_{p}!}$)

  * •

$\mathcal{N}_{0,\sigma}$ accounts for the read-out noise.

While this model seemingly involves many parameters ($G_{2},\ \lambda,\
\sigma$), all of these can actually be expressed in terms of $G_{2}$ through
the following equations

$\begin{split}\frac{\lambda G_{2}}{4}&\simeq\text{mean of empirical data }\\\
\frac{2G_{2}^{2}\lambda}{4^{2}}+\sigma^{2}&\simeq\text{ variance of empirical
data }\\\ \end{split}$ (15)

In practice, the mean gain over all experimental series is
$\overline{G_{2}}=23.46$. The standard deviation of the mean gain
$\overline{G_{2}}$ for all experimental series is 0.48. The standard deviation
of $G_{2}$ within a single experimental series is $0.8$. The typical fit
standard error is $0.22$. These values are in very good agreement with the
results of the first method.

For this value of $G_{2}$, and our camera having 590 amplification stages,
validity of eq.(12), hence of eq.(13), is expected to be better than $99\%$
according to formula (8) in Robbins and Hadwen (2003).

Figure 5: Black points : experimental distribution of the number of counts per
pixel. Blue line : Fit of the data. This single parameter fit directly gives
the value of $G_{2}=24.02$ for the image considered here

Finally, we give the expression of the contribution of the shot noise on the
standard deviation of normalized magnetization $M_{z}$ in our system. With
$S_{m_{s}}$ the signal for the spin components $m_{s}$,
$M_{z}=\sum_{m_{s}}m_{s}S_{m_{S}}/\sum_{m_{s}}S_{m_{S}}$, one gets, assuming
independent noises:

$\sigma_{SN}=\frac{\sqrt{2G_{exp}\sum_{m_{s}}m_{s}^{2}S_{m_{S}}}}{\sum_{m_{s}}S_{m_{S}}}$
(16)

The corresponding variance is proportional to the total number of atoms $N$,
and inversely proportional to the photon collection efficiency.

The shot noise has a significant impact on the observed magnetization
fluctuations, as described in the main paper. In our analysis, we thus
subtract from the experimental variance the sum of the estimated variance due
to shot noise and of the estimated variance due to fit uncertainty. The latter
is then estimated assuming that the noise on each pixel of the image is
independent of the signal at this pixel.

As shot noise can be non-negligible and is signal-dependent, it is not a
priori justified to assume a signal-independent noise to deduce the fit-noise.
We therefore have also estimated the combined effect of the shot noise and
other signal-independent noises (such as read noise) in the following manner:
we estimate the variance of the fit, using a $\chi^{2}$ that is now normalized
to the noise on each pixel. This latter noise is obtained by adding in
quadrature the signal-independent noise (obtained by statistics on the parts
of the camera where the signal is negligible) and the estimated shot noise per
pixel (see above). We find that both methods give very similar results, which
validate the approach that we follow in the main part of the paper.

###  I.5 Details of numerical simulations and determination of the quadratic
field value

Here we provide details regarding our numerical simulations, as well as the
determination of the quadratic light shift $B_{Q}$ induced by the lattice
lasers, which is the only free parameter used in our numerical simulations in
the main text. To examine the population and correlation dynamics for large
system sizes, we make use of the generalized discrete truncated Wigner
approximation (GDTWA), previously introduced in Zhu  _et al._ (2019).

In order to make a comparison with experimental results, we must first
determine the best-fit value of the quadratic light shift, which is known to
be present in the experiment. We perform numerical GDTWA simulations of the
Hamiltonian in the main text with an added term of the form
$B_{Q}\sum_{i}(\hat{s}_{z}^{i})^{2}$. For each value of $B_{Q}$ that we
consider, we compute the quantity
$\chi^{2}_{m_{s}}(t)=(p_{m_{s}}^{\left[\mathit{sim}\right]}(t)-p_{m_{s}}^{\left[\mathit{exp}\right]}(t))^{2}/\sigma^{2}_{m_{s}}(t)$,
where $p_{m_{s}}^{\left[\mathit{sim/exp}\right]}(t)$ denotes the
simulated/experimental values of $p_{m_{s}}(t)$, and $\sigma^{2}_{m_{s}}(t)$
corresponds to the experimental error (we assume comparatively negligible
GDTWA sampling error). We then average over all $m_{s}$ and available $t$ to
obtain the mean-squared error $\chi^{2}=\overline{\chi^{2}_{m_{s}}(t)}$.
Minimizing this quantity over $B_{Q}$, we find an optimal value of
$B_{Q}\approx-5.1$ Hz with $\chi^{2}=15.3$. In Fig. 6, we plot $\chi^{2}$ over
the range of $B_{Q}$ we consider, and compare the resulting dynamics of
$p_{m_{s}}$ for the best-fit $B_{Q}$. We observe that our simulations provide
decent agreement with the experimental results, capturing the relaxation
timescales and steady-state values of $p_{m_{s}}$.

All of our results are obtained for a lattice size of $L_{x}\times L_{y}\times
L_{z}=5\times 3\times 5$ with periodic boundary conditions, which we have
verified is enough to produce results for $p_{m_{s}}(t)$ that are convergent
in system size along each dimension, within relevant timescales and
experimental error bars.

Figure 6: Determination of the quadratic light shift from the experimental
spin populations. (Left) Value of $\chi^{2}=\overline{\chi_{m_{s}}^{2}}$
(orange line) for various $B_{Q}$. We also plot $\chi^{2}_{m_{s}}$ for each
$m_{s}$ averaged over all available times; dashed (solid) lines denote
positive (negative) $m_{s}$. Vertical dotted line denotes optimal value of
$B_{Q}$ (minimum of $\chi^{2}$). (Right) GDTWA results for best-fit
$B_{Q}\approx-5.1$ Hz, compared with experimental results for the population
dynamics; positive (negative) $m_{s}$ displayed in lighter (darker) shading.

###  I.6 Estimate of the correlation length from the measurement of $C_{z}$
and $\Sigma_{z}$

In this section we provide an estimate of the characteristic correlation
length associated with the dynamical onset of correlations in the system,
starting from the measurement of $C_{z}$ and $\Sigma_{z}$.

To this scope, let us introduce the spin-spin correlation function
$G_{z}(i,j)=\langle s_{z}^{i}s_{z}^{j}\rangle-\langle s_{z}^{i}\rangle\langle
s_{z}^{j}\rangle$. When focusing on the long-time regime, in which the system
is expected to thermalize, we can in general expect $G_{z}$ to be short-
ranged, namely decaying exponentially with the distance between the sites,
with a decay rate given by the correlation length $\xi$. The spatial
anisotropy of the dipolar interactions, as well as of the optical lattice used
in the experiment, would suggest that there are in fact several correlation
lengths when moving along different lattice directions; nonetheless for
simplicity we shall neglect this aspect, and assume that correlations are
spatially isotropic.

Moreover, in the same spirit we shall assume that the system is
translationally invariant, namely that $G_{z}(i,j)=G_{z}(r_{ij})$, where
$r_{ij}$ is the distance between the $i$-th and $j$-th spin. $G_{z}(r)$ can
therefore be taken to be $|G(z)|\approx G_{z}(0)\exp(-r/\xi)$. The zero-range
correlations $G_{z}(0)$ coincide with the on-site spin fluctuations, namely
with $\Sigma_{z}/N$. On the other hand $C_{z}$ is the sum of the offsite
correlations; given that $C_{z}$ is negative, it is plausible to assume that
the sign of all the correlations for $r>0$ is globally negative. Putting all
these assumptions together, we posit the simple functional form
$G_{z}(r>0)\approx-G_{z}(0)e^{-r/\xi}$ for the off-site correlations.

By definition of $C_{z}$ we have that

$\frac{C_{z}}{N}=\sum_{\bm{r}\neq
0}G_{z}(|\bm{r}|)\approx\frac{\Sigma_{z}}{N}\sum_{\bm{r}\neq 0}e^{-r/\xi}$
(17)

from which we deduce the integral relationship between $C_{z}$, $\Sigma_{z}$
and $\xi$

$\sum_{\bm{r}\neq 0}e^{-r/\xi}=\frac{|C_{z}|}{\Sigma_{z}}~{}.$ (18)

At long times ($t=100$ ms) we measure $C_{z}/N\approx-1.3$ and
$\Sigma_{z}/N\approx 2.8$. Solving Eq. (18) for $\xi$ numerically, we obtain
the value $\xi=0.315$ quoted in the main text.

Let us remark that the working assumption $G_{z}(r>0)\approx-
G_{z}(0)e^{-r/\xi}$, which attributes the same amplitude to the on-site
fluctuations and to the exponential tail of the correlations, is actually
assuming a maximum concentration of correlations at short distances. Therefore
its use provides in practice a _lower bound_ to the correlation length $\xi$.
A more general form would have been $G_{z}(r>0)\approx-Ae^{-r/\xi}$ with
$A\leq G_{z}(0)$, leading to a larger $\xi$ estimate. Yet the experimental
data at hand do not allow us to estimate $\xi$ and $A$ independently.

###  I.7  Effect of losses on correlations

To estimate the effect of losses on the measured covariances and correlations,
we have developed a simple heuristic statistical model. This model ignores the
effect of the dipole-dipole interactions between atoms that is described by
the secular Hamiltonian in the paper, but it includes dipolar relaxation. It
considers the case of a system where all spin are initially tilted by
$\theta=\pi/2$ compared to the magnetic field.

Consider the initial number of single- and double-occupied sites, denoted $S$
and $D$, respectively. The probability that a given atom initially occupies a
single-occupied site is then $P_{S}=S/(S+2D)$, while the probability that it
occupies a doubly occupied state is $P_{D}=2D/(S+2D)$. Since the magnetic
field is sufficiently large, dipolar relaxation only affects doubly-occupied
sites; for an atom occupying such a site and with a spin projection $m_{i}$,
the inverse lifetime due to dipolar relaxation is
$\Gamma_{m_{i}}=\sum_{m_{s}}p_{m_{s}}\Gamma_{m_{i},m_{m_{s}}}$. Here,
$p_{m_{s}}$ are the fractional populations, which in principle are time-
dependent but we will take to equal the populations in the initial state,
since dipolar relaxation lead to small dynamics in the fractional populations
Kechadi (2019). The $\Gamma_{m_{i},m_{m_{s}}}$ are calculated in the Born
approximation Pasquiou  _et al._ (2010); Kechadi (2019), and are spin-
dependent through the angular terms in the spin operators. Thus, the
probability to measure a given atom in the state $m_{i}$ can be written as

$\displaystyle
P_{m_{i}}=\alpha\frac{p_{m_{i}}}{S+2D}\times\left[S+2De^{-\Gamma_{m_{i}}t}\right].$

The time-dependent constant $\alpha$ is defined by the normalization condition
$\sum_{m_{i}}P_{m_{i}}=1$.

In the absence of interactions, two atoms on separate lattice sites will decay
independently of each other. Since the probability that any two given atoms
will occupy the same lattice site is very small, we make the additional
assumption that _any_ two atoms will decay independently of each other. The
measurement of the spin projection $m_{i}$ over $N$ uncorrelated atoms is then
described by a multinomial process with $N$ independent trials and parameters
$P_{m_{i}}$ describing the probability of measurement outcome $m_{i}$ for each
atom. Thus, the variance and covariance are given by

$\displaystyle{\rm Var}[N_{m_{i}}]=NP_{m_{i}}(1-P_{m_{i}})$ (19)
$\displaystyle{\rm Cov}[N_{m_{i}},N_{m_{j}}]=-NP_{m_{i}}P_{m_{j}}.$ (20)

This allows for an estimate of ${\rm Var}[S_{z}]$. As expected we find that at
$t=0$ (losses have not occurred) and at large $t$ (where all doublons have
disappeared), the outcome is identical, and equal to the expected variance
$3/2N$. For intermediate times, the variance varies by typically less than 5
percent. This indicates that the effect of losses on correlations may safely
be neglected in our experiment.

###  I.8 Study of the stationarity of the preparation noise through the whole
dynamics

In the case that there are fluctuations in the initial rotation angle
$\theta$, the total variances can be obtained as

$\displaystyle{\rm Var}^{\rm tot}[O]$ $\displaystyle=$ $\displaystyle\int
d\theta p(\theta)\langle\hat{O}^{2}\rangle_{\theta}-\int d\theta
p(\theta)\langle\hat{O}\rangle_{\theta}\int d\theta
p(\theta)\langle\hat{O}\rangle_{\theta},$

where $p(\theta)$ is the distribution of the initial rotation angles in the
presence of rf noise.

It is straightforward to see the above is equivalent to

$\displaystyle{\rm Var}^{\rm tot}[O]$ $\displaystyle=$
$\displaystyle\overline{{\rm
Var}[O]_{\theta}}+\overline{(\langle\hat{O}\rangle_{\theta}-\overline{\langle\hat{O}\rangle_{\theta}})^{2}}$
(22) $\displaystyle=\int d\theta
p(\theta)[\langle\hat{O}^{2}\rangle_{\theta}-\langle\hat{O}\rangle^{2}_{\theta}]$
$\displaystyle+$ $\displaystyle[\int d\theta
p(\theta)\langle\hat{O}\rangle_{\theta}^{2}-(\int d\theta
p(\theta)\langle\hat{O}\rangle_{\theta})^{2}],$

where $\overline{O}$ denotes the average over the distribution of $\theta$,
and $\langle\hat{O}\rangle$ denotes the quantum expectation value. That is, we
can compute the total variance either using Eq. (LABEL:eq:totvar1) or directly
using Eq. (22). In the following, we will further show some simplified results
from Eq. (22).

In the case that $p(\theta)$ is a narrow distribution around a certain angle
$\theta_{0}$, eg.

$\displaystyle p(\theta)$ $\displaystyle\propto$ $\displaystyle
e^{-\frac{(\theta-\theta_{0})^{2}}{2\sigma^{2}_{\theta}}},$ (23)

with $\sigma_{\theta}\ll 1$, we can expand the functions inside the integral
$\int d\theta p(\theta)F(\theta)$ as

$\displaystyle F(\theta)\\!$ $\displaystyle\approx$
$\displaystyle\\!F(\theta_{0})+\\!\frac{\partial
F}{\partial\theta}\bigg{|}_{\theta_{0}}\\!\\!\\!(\theta-\theta_{0})+\frac{\partial^{2}F}{\partial\theta^{2}}\bigg{|}_{\theta_{0}}\\!\\!\\!\frac{(\theta-\theta_{0})^{2}}{2}\\!+\\!\mathcal{O}(\theta-\theta_{0})^{3}.$

Substitute this into either of the above equations for the total
(co)variances, one obtains

$\displaystyle{\rm Var}^{\rm tot}[O]$ $\displaystyle\approx$
$\displaystyle\int d\theta p(\theta){\rm
Var}[O]_{\theta}+\left[\frac{\partial\langle\hat{O}\rangle}{\partial\theta}\bigg{|}_{\theta_{0}}\right]^{2}\sigma_{\theta}^{2}$
$\displaystyle\approx$ $\displaystyle{\rm
Var}[O]_{\theta_{0}}+\frac{1}{2}\frac{\partial^{2}{\rm
Var}[O]}{\partial\theta^{2}}\bigg{|}_{\theta_{0}}\sigma_{\theta}^{2}+\left[\frac{\partial\langle\hat{O}\rangle}{\partial\theta}\bigg{|}_{\theta_{0}}\right]^{2}\sigma_{\theta}^{2},$

where ${\rm Var}[O]$ is the quantum variance from an initial state $\theta$,
as calculated in the previous section. Note, for the initial state considered
in experiment, the first term ${\rm Var}[O]_{\theta_{0}}$ is $\propto N$, and
the second term $\propto N\sigma_{\theta}^{2}$, while the last term in Eq.
(I.8) is $\propto N^{2}\sigma_{\theta}^{2}$. So when $N\gg 1$, the second term
is negligible compared to the last term, and can be dropped. This means
roughly one can estimate the total variance as

$\displaystyle{\rm Var}^{\rm tot}[O]$ $\displaystyle\approx$
$\displaystyle{\rm
Var}[O]_{\theta_{0}}+\left[\frac{\partial\langle\hat{O}\rangle}{\partial\theta}\bigg{|}_{\theta_{0}}\right]^{2}\sigma_{\theta}^{2},$
(27)

ie., the simple summation of the contribution from two quadratures. To have
quantum noise that is significant compared to the technical noise
contribution, this also suggests that a small rf noise is needed
$\sigma_{\theta}^{2}\sim 1/N$.

###  I.9 Comparison between intrasite and intersite spin correlations

Figure 7: Spin correlations for doublons, $\Sigma_{aa}$ (blue), $\Sigma_{ab}$
(orange), and $\mathcal{C}_{z}$ (green). Results obtained from GDTWA
simulations for a 3D lattice with $N_{d}=1372$, with $B_{Q}=-5$Hz and zero
magnetic field gradients. The lattice sites are all doubly occupied, with
total atom number $N=2N_{d}$.

As discussed in Sec. I.2, in the presence of doublons, in addition to the
intra-spin correlations $\Sigma_{z}$, correlations can build up between atoms
in the same site,
$\Sigma_{ab}=2\sum_{i}\langle\hat{s}_{z}^{i1}\hat{s}_{z}^{i2}\rangle$, as well
as between different sites, $\mathcal{C}_{z}=\sum_{i,j\neq
i}^{\rm\\#of~{}sites}\langle\hat{s}_{z}^{i}\hat{s}_{z}^{j}\rangle$. Describing
doublons as $S=6$ particles, we find at short time, the growth of these
correlations takes the form

$\displaystyle\Sigma_{ab}$ $\displaystyle=$
$\displaystyle-27t^{2}(4B_{Q}\overline{V}-3V_{\rm eff})^{2},$ (28)
$\displaystyle\mathcal{C}_{il}$ $\displaystyle=$
$\displaystyle\frac{99t^{2}}{2}(4B_{Q}\overline{V}-3V_{\rm eff}^{2}),$ (29)

where $\overline{V}=\sum_{j<k}V_{jk}/2N_{d}$ and $V_{\rm
eff}^{2}=\sum_{j<k}V_{jk}^{2}/2N_{d}$, with the summation running over
different lattice sites populated by doublons. That is, the intra-site and
inter-site correlations have different growth rates and signs. As shown in
Sec. I.2, $\Sigma_{aa}$ can be obtained from experiment with collective
measurement. We have verified with numerical simulations that it can be
related to $\Sigma_{ab}$ via

$\displaystyle\Sigma_{ab}\\!=\\!\frac{6}{5}\left(\Sigma_{aa}(t)-\frac{2}{7+5\cos(2\theta)}\Sigma_{aa}(0)\right),$
(30)

where $\theta$ is the tipping angle of the initial state, and in this work
$\theta=\pi/2$. Then $\mathcal{C}_{z}$ can be obtained from
$\mathcal{C}_{z}=3N/2-\Sigma_{aa}-\Sigma_{ab}$.

In Fig. 7 we use GDTWA to calculate the dynamics of different correlations up
to a long timescale relevant for experiment, which shows that atoms in
different sites become significantly anti-correlated under the dipolar
interactions in our system.

## References

  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Cavalcanti and Skrzypczyk (2016) D. Cavalcanti and P. Skrzypczyk, Reports on Progress in Physics 80, 024001 (2016).
  * Pezzè  _et al._ (2018) L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Rev. Mod. Phys. 90, 035005 (2018).
  * Brunner  _et al._ (2014) N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Rev. Mod. Phys. 86, 419 (2014).
  * Ollivier and Zurek (2001) H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).
  * Luo (2008) S. Luo, Phys. Rev. A 77, 042303 (2008).
  * Bera  _et al._ (2017) A. Bera, T. Das, D. Sadhukhan, S. S. Roy, A. Sen(De), and U. Sen, Reports on Progress in Physics 81, 024001 (2017).
  * Wang  _et al._ (2020) Y. Wang, Z. Hu, B. C. Sanders, and S. Kais, Frontiers in Physics 8, 479 (2020).
  * Ferraro  _et al._ (2010) A. Ferraro, L. Aolita, D. Cavalcanti, F. M. Cucchietti, and A. Acín, Phys. Rev. A 81, 052318 (2010).
  * Kaufman  _et al._ (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * Fadel  _et al._ (2018) M. Fadel, T. Zibold, B. Décamps, and P. Treutlein, Science 360, 409 (2018).
  * Kunkel  _et al._ (2018) P. Kunkel, M. Prüfer, H. Strobel, D. Linnemann, A. Frölian, T. Gasenzer, M. Gärttner, and M. K. Oberthaler, Science 360, 413 (2018).
  * Lange  _et al._ (2018) K. Lange, J. Peise, B. Lücke, I. Kruse, G. Vitagliano, I. Apellaniz, M. Kleinmann, G. Tóth, and C. Klempt, Science 360, 416 (2018).
  * Schmied  _et al._ (2016) R. Schmied, J.-D. Bancal, B. Allard, M. Fadel, V. Scarani, P. Treutlein, and N. Sangouard, Science 352, 441 (2016).
  * Sørensen and Mølmer (2001) A. S. Sørensen and K. Mølmer, Phys. Rev. Lett. 86, 4431 (2001).
  * Vitagliano  _et al._ (2011) G. Vitagliano, P. Hyllus, I. n. L. Egusquiza, and G. Tóth, Phys. Rev. Lett. 107, 240502 (2011).
  * Vitagliano  _et al._ (2014) G. Vitagliano, I. Apellaniz, I. n. L. Egusquiza, and G. Tóth, Phys. Rev. A 89, 032307 (2014).
  * Bookjans  _et al._ (2011) E. M. Bookjans, C. D. Hamley, and M. S. Chapman, Phys. Rev. Lett. 107, 210406 (2011).
  * Gross  _et al._ (2011) C. Gross, H. Strobel, E. Nicklas, T. Zibold, N. Bar-Gill, G. Kurizki, and M. K. Oberthaler, Nature 480, 219 (2011).
  * Lücke  _et al._ (2011) B. Lücke, M. Scherer, J. Kruse, L. Pezzé, F. Deuretzbacher, P. Hyllus, O. Topic, J. Peise, W. Ertmer, J. Arlt, L. Santos, A. Smerzi, and C. Klempt, Science 334, 773 (2011).
  * Qu  _et al._ (2020) A. Qu, B. Evrard, J. Dalibard, and F. Gerbier, Phys. Rev. Lett. 125, 033401 (2020).
  * Linnemann  _et al._ (2016) D. Linnemann, H. Strobel, W. Muessel, J. Schulz, R. J. Lewis-Swan, K. V. Kheruntsyan, and M. K. Oberthaler, Phys. Rev. Lett. 117, 013001 (2016).
  * Evrard  _et al._ (2021) B. Evrard, A. Qu, J. Dalibard, and F. Gerbier, Science 373, 1340 (2021).
  * Law  _et al._ (1998) C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett. 81, 5257 (1998).
  * Lepoutre  _et al._ (2019) S. Lepoutre, J. Schachenmayer, L. Gabardos, B. H. Zhu, B. Naylor, E. Maréchal, O. Gorceix, A. M. Rey, L. Vernac, and B. Laburthe-Tolra, Nature Communications 10, 1714 (2019).
  * Patscheider  _et al._ (2020) A. Patscheider, B. Zhu, L. Chomaz, D. Petter, S. Baier, A.-M. Rey, F. Ferlaino, and M. J. Mark, Phys. Rev. Research 2, 023050 (2020).
  * (27) See Supplemental Material.
  * Fersterer  _et al._ (2019) P. Fersterer, A. Safavi-Naini, B. Zhu, L. Gabardos, S. Lepoutre, L. Vernac, B. Laburthe-Tolra, P. B. Blakie, and A. M. Rey, Phys. Rev. A 100, 033609 (2019).
  * Lepoutre  _et al._ (2018a) S. Lepoutre, K. Kechadi, B. Naylor, B. Zhu, L. Gabardos, L. Isaev, P. Pedri, A. M. Rey, L. Vernac, and B. Laburthe-Tolra, Phys. Rev. A 97, 023610 (2018a).
  * Lepoutre  _et al._ (2018b) S. Lepoutre, L. Gabardos, K. Kechadi, P. Pedri, O. Gorceix, E. Maréchal, L. Vernac, and B. Laburthe-Tolra, Phys. Rev. Lett. 121, 013201 (2018b).
  * de Paz  _et al._ (2013) A. de Paz, A. Chotia, E. Maréchal, P. Pedri, L. Vernac, O. Gorceix, and B. Laburthe-Tolra, Phys. Rev. A 87, 051609(R) (2013).
  * Zhu  _et al._ (2019) B. Zhu, A. M. Rey, and J. Schachenmayer,  , 45 (2019), arXiv:1905.08782 .
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016).
  * (34) S. M. Bihui Zhu and A. M. Rey, Work in progress .
  * Ammann and Christensen (1997) H. Ammann and N. Christensen, Phys. Rev. Lett. 78, 2088 (1997).
  * Robbins and Hadwen (2003) M. Robbins and B. Hadwen, Electron Devices, IEEE Transactions on 50, 1227 (2003).
  * Kechadi (2019) K. Kechadi, PhD Thesis (2019).
  * Pasquiou  _et al._ (2010) B. Pasquiou, G. Bismut, Q. Beaufils, A. Crubellier, E. Maréchal, P. Pedri, L. Vernac, O. Gorceix, and B. Laburthe-Tolra, Phys. Rev. A 81, 042716 (2010).


# Mixing ADAM and SGD: a Combined Optimization Method

Nicola Landro, Ignazio Gallo, Riccardo La Grassa
University of Insubria
{nlandro, ignazio.gallo<EMAIL_ADDRESS>
###### Abstract

Optimization methods (optimizers) get special attention for the efficient
training of neural networks in the field of deep learning. In literature there
are many papers that compare neural models trained with the use of different
optimizers. Each paper demonstrates that for a particular problem an optimizer
is better than the others but as the problem changes this type of result is no
longer valid and we have to start from scratch. In our paper we propose to use
the combination of two very different optimizers but when used simultaneously
they can overcome the performances of the single optimizers in very different
problems. We propose a new optimizer called MAS (Mixing ADAM and SGD) that
integrates SGD and ADAM simultaneously by weighing the contributions of both
through the assignment of constant weights. Rather than trying to improve SGD
or ADAM we exploit both at the same time by taking the best of both. We have
conducted several experiments on images and text document classification,
using various CNNs, and we demonstrated by experiments that the proposed MAS
optimizer produces better performance than the single SGD or ADAM optimizers.
The source code and all the results of the experiments are available online at
the following link https://gitlab.com/nicolalandro/multi_optimizer

##  1 Introduction

Stochastic Gradient Descent [1] (SGD) is the dominant method for solving
optimization problems. SGD iteratively updates the model parameters by moving
them in the direction of the negative gradient calculated on a mini-batch
scaled by the step length, typically referred to as the learning rate. It is
necessary to decay this learning rate as the algorithm proceeds to ensure
convergence. Manually adjusting the learning rate decay in SGD is not easy. To
address this problem, several methods have been proposed that automatically
reduce the learning rate. The basic intuition behind these approaches is to
adaptively tune the learning rate based only on recent gradients; therefore,
limiting the reliance on the update to only a few past gradients. ADAptive
Moment estimation [2] (ADAM) is one of several methods based on this update
mechanism [3]. On the other hand, adaptive optimization methods such as ADAM,
even though they have been proposed to achieve a rapid training process, are
observed to generalize poorly with respect to SGD or even fail to converge due
to unstable and extreme learning rates [4]. To try to overcome the problems of
both of these types of optimizers and at the same time try to exploit their
advantages, we propose an optimizer that combines them in a new optimizer.

Figure 1: Intuitive representation of the idea behind the proposed MAS (Mixing
ADAM and SGD) optimizer: the weights are modified simultaneously by both the
optimizers.

As depicted in Figure 1, the basic idea of the MAS optimizer here proposed, is
to combine two different known optimizers and automatically go quickly towards
the direction of both on the surface of the loss function when the two
optimizers agree (see geometric example in Figure 2a). While when the two
optimizers used in the combination do not agree, our solution always goes
towards the predominant direction between the two but slowing down the speed
(see example of Figure 2b).

Analyzing the literature there are many papers that compare neural models
trained with the use of different optimizers [5, 6, 7, 8] or that propose
modifications for existing optimizers [9, 10, 11], always aimed at improving
the results on a subset of problems. Each paper demonstrates that an optimizer
is better than the others but as the problem changes this type of result is no
longer valid and we have to start from scratch. In our paper we propose to
combine simultaneously two different optimizers like SGD and ADAM, to overcome
the performances of the single optimizers in very different problems.

Contributions. In light of the above, we set out the main contributions of our
paper.

  * •

We demonstrate experimentally that the combination of different optimizers in
a new optimizer that incorporates them leads to a better generalization
capacity in different contexts.

  * •

We apply the proposed solution to well-known computer vision and text analysis
problems and in both these domains, we obtain excellent results, demonstrating
that our solution exceeds the generalization capacity of the starting ADAM and
SGD optimizers.

  * •

We open the way to this new type of optimizers with which it will be possible
to exploit the positive aspects of many existing optimizers, combining them
with each others to build new and more efficient optimizers.

  * •

To facilitate the understanding of the MAS optimizer and to allow other
researchers to be able to run the experiments and extend this idea, we release
the source code and setups of the experiments at the following URL [12]

##  2 Related Work

In the literature, there aren’t many papers that try to combine different
optimizers together. In this section, we report some of the more recent papers
that in some ways use different optimizers in the same learning process.

In [13] the authors investigate a hybrid strategy, called SWATS, which starts
training with an adaptive optimization method and switches to SGD when
appropriate. This idea starts from the observation that despite superior
training results, adaptive optimization methods such as ADAM generalize poorly
compared to SGD because they tend to work well in the early part of the
training, but are overtaken by SGD in the later stages of training. In
concrete terms, SWATS is a simple strategy that goes from Adam to SGD when an
activation condition is met. The experimental results obtained in this paper
are not so different from ADAM or SGD used individually, so the authors
concluded that using SGD with perfect parameters is the best idea. In our
proposal, we want to combine two well-known optimizers to create a new one
that uses simultaneously two different optimizers from the beginning to the
end of the training process.

ESGD is a population-based Evolutionary Stochastic Gradient Descent framework
for optimizing deep neural networks [14]. In this approach, individuals in the
population optimized with various SGD-based optimizers using distinct hyper-
parameters are considered competing species in a context of coevolution. The
authors experimented with optimizer pools consisting of SGD and ADAM variants
where it is often observed that ADAM tends to be aggressive early on, but
stabilizes quickly, while SGD starts slowly but can reach a better local
minimum. ESGD can automatically choose the appropriate optimizers and their
hyper-parameters based on the fitness value during the evolution process so
that the merits of SGD and ADAM can be combined to seek a better local optimal
solution to the problem of interest. In the method we propose, we do not need
another approach, such as the evolutionary one, to decide which optimizer to
use and with which hyper-parameters, but it is the same approach that decides
at each step what is the contribution of SGD and that of ADAM.

In this paper, we also compare our MAS optimizer with ADAMW [15] [16], which
is a version of ADAM in which weight decay is decoupled from $L_{2}$
regularization. This optimizer offers good generalization performance,
especially for text analysis, and since we also perform some experimental
tests on text classification, then we also compare our optimizer with ADAMW.
In fact, ADAMW is often used with BERT [17] applied to well-known datasets for
text classification.

Figure 2: Graphical representation of the basic idea for the proposed MAS
optimizer. In (a), if the two translations $\vec{w}_{1}$ and $\vec{w}_{2}$
obtained from two different optimizers are similar, then the resulting
translation $\vec{w}_{1}+\vec{w}_{2}$ is boosted. In (b), if the translations
$\vec{w}_{1}$ and $\vec{w}_{2}$ go in two different directions, then the
resulting translation is smaller. We also use two hyper-parameters
$\lambda_{1}$ and $\lambda_{2}$ to weigh the contribution of the two
optimizers. Algorithm 1 Stochastic Gradient Descent (SGD)

Input: the weights $w_{k}$, learing rate $\eta$, weight decay $\gamma$,
dampening $d$, boolean $nesterov$

1:$v_{0}=0$

2:function $\Delta_{\text{SGD}}$($w_{k}$, $\nabla$, $\gamma$, $\mu$, $d$,
$nesterov$)

3: $\widehat{\nabla}=\nabla+w_{k}\cdot\gamma$

4: if $m\neq 0$ then

5: if $k=0$ then

6: $v_{k}=\widehat{\nabla}$

7: else

8: $v_{k}=v_{k-1}\cdot\mu+\widehat{\nabla}\cdot(1-d)$

9: end if

10: if $nesterov=True$ then

11: $v_{k}=\widehat{\nabla}+v_{k}\cdot\mu$

12: end if

13: end if

14: return $v_{k}$

15:end function

16:for batches do

17:
$w_{k+1}=w_{k}-\eta_{s}\cdot\Delta_{\text{SGD}}(w_{k},\nabla,\gamma,\mu,d,nesterov)$

18:end for

##  3 Preliminaries

Training neural networks is equivalent to solving the following optimization
problem:

$\min_{w\in\mathbb{R}^{n}}\mathcal{L}(w)$ (1)

where $\mathcal{L}$ is a loss function and $w$ are the weights. The iterations
of an SGD [1] optimizer can be described as:

$w_{k+1}=w_{k}-\eta\cdot\nabla\mathcal{L}(w)$ (2)

where $w_{k}$ denotes the weights $w$ at the $k$-th iteration, $\eta$ denote
the learning rate and $\nabla\mathcal{L}(w)$ denotes the stochastic gradient
calculated at $w_{k}$. To propose a stochastic gradient calculated as generic
as possible, we introduce the weight decay [18] strategy, often used in many
SGD implementations. The weight decay can be seen as a modification of the
$\nabla\mathcal{L}(w)$ gradient, and in particular, we can describe it as
follows:

$\widehat{\nabla}\mathcal{L}(w_{k})=\nabla\mathcal{L}(w_{k})+w_{k}\cdot\gamma$
(3)

where $\gamma$ is a small scalar called weight decay. We can observe that if
the weight decay $\gamma$ is equal to zero then
$\widehat{\nabla}\mathcal{L}(w)=\nabla\mathcal{L}(w)$. Based on the above, we
can generalize the Eq. 2 to the following one that include weight decay:

$w_{k+1}=w_{k}-\eta\cdot\widehat{\nabla}\mathcal{L}(w)$ (4)

The SGD algorithm described up to here is usually used in combination with
momentum, and in this case, we refer to it as SGD(M) [19] (Stochastic Gradient
Descend with Momentum). SGD(M) almost always works better and faster than SGD
because the momentum helps accelerate the gradient vectors in the right
direction, thus leading to faster convergence. The iterations of SGD(M) can be
described as follows:

$v_{k}=\mu\cdot v_{k-1}+\widehat{\nabla}\mathcal{L}(w)$ (5)
$w_{k+1}=w_{k}-\eta\cdot v_{k}$ (6)

where $\mu\in[0,1)$ is the momentum parameter and for $k=0$, $v_{0}$ is
initialized to $0$. The simpler methods of momentum have an associated damping
coefficient [20], which controls the rate at which the momentum vector decays.
The dampening coefficient changes the momentum as follow:

$v_{d_{k}}=m\cdot v_{k-1}+\widehat{\nabla}\mathcal{L}(w)\cdot(1-d)$ (7)

where $0\leq d<1$, is the dampening value, so the final SGD with momentum and
dampening coefficients can be seen as follow:

$w_{k+1}=w_{k}-\eta\cdot v_{d_{k}}$ (8)

Nesterov momentum [21] is an extension of the moment method that approximates
the future position of the parameters that takes into account the movement.
The SGD with nesterov transforms again the $v_{k}$ of Eq. 5, more precisely:

$v_{n_{k}}=\widehat{\nabla}\mathcal{L}(w)+v_{d_{k}}\cdot m$ (9)
$w_{k+1}=w_{k}-\eta\cdot v_{n_{k}}$ (10)

The complete SGD algorithm, used in this paper, is showed in Alg. 1.

Algorithm 2 ADaptive Moment Estimation (ADAM)

Input: the weights $w_{k}$, learing rate $\eta$, weight decay $\gamma$,
$\beta_{1}$, $\beta_{2}$, $\epsilon$, boolean $amsgrad$

1:$m_{0}=0$

2:$v^{a}_{0}=0$

3:$\widehat{v}_{0}=0$

4:function $\Delta_{\text{ADAM}}$($w_{k}$, $\nabla$, $\eta$, $\gamma$,
$\beta_{1}$, $\beta_{2}$, $\epsilon$, $amsgrad$)

5: $\widehat{\nabla}=\nabla+w_{k}\cdot\gamma$

6: $m_{k}=m_{k-1}\cdot\beta_{1}+\widehat{\nabla}\cdot(1-\beta_{1})$

7:
$v^{a}_{k}=v^{a}_{k-1}\cdot\beta_{2}+\widehat{\nabla}\cdot\widehat{\nabla}\cdot(1-\beta_{2})$

8: if $amsgrad=True$ then

9: $\widehat{v}_{k}=\text{max}(\widehat{v}_{k-1},v^{a}_{k})$

10: $denom=\frac{\sqrt{\widehat{v}_{k}}}{\sqrt{1-\beta_{2}}+\epsilon}$

11: else

12: $denom=\frac{\sqrt{v_{k}^{a}}}{\sqrt{1-\beta_{2}}+\epsilon}$

13: end if

14: $\eta_{a}=\frac{\eta}{1-\beta_{1}}$

15: $d_{k}=\frac{m_{k}}{denom}$

16: return $d_{k},\eta_{a}$

17:end function

18:for batches do

19:
$d_{k},\eta_{a}=\Delta_{\text{ADAM}}(w_{k},\nabla,\eta,\gamma,\beta_{1},\beta_{2},\epsilon,amsgrad)$

20: $w_{k+1}=w_{k}-\eta_{a}\cdot d_{k}$

21:end for

ADAM [2] (ADAptive Moment estimation) optimization algorithm is an extension
to SGD that has recently seen broader adoption for deep learning applications
in computer vision and natural language processing. ADAM’s equation for
updating the weights of a neural network by iterating over the training data
can be represented as follows:

$m_{k}=\beta_{1}\cdot
m_{k-1}+(1-\beta_{1})\cdot\widehat{\nabla}\mathcal{L}(w_{k})$ (11)
$v_{k}^{a}=\beta_{2}\cdot
v_{k-1}+(1-\beta_{2})\cdot\widehat{\nabla}\mathcal{L}(w_{k})^{2}$ (12)
$w_{k+1}=w_{k}-\eta\cdot\frac{\sqrt{1-\beta_{2}}}{1-\beta_{1}}\cdot\frac{m_{k}}{\sqrt{v_{k}^{a}}+\epsilon}$
(13)

where $m_{k}$ and $v_{k}^{a}$ are estimates of the first moment (the mean) and
the second moment (the non-centered variance) of the gradients respectively,
hence the name of the method. $\beta_{1}$, $\beta_{2}$ and $\epsilon$ are
three new introduced hyper-parameters of the algorithm. AMSGrad [22] is a
stochastic optimization method that seeks to fix a convergence issue with Adam
based optimizers. AMSGrad uses the maximum of past squared gradients $v_{k-1}$
rather than the exponential average to update the parameters:

$\widehat{v}_{k}=\text{max}(\widehat{v}_{k-1},v_{k}^{a})$ (14)

$w_{k+1}=w_{k}-\eta\cdot\frac{\sqrt{1-\beta_{2}}}{1-\beta_{1}}\cdot\frac{m_{k}}{\sqrt{\widehat{v}_{k}}+\epsilon}$
(15)

The complete ADAM algorithm, used in this paper, is showed in Alg. 2.

##  4 Proposed Approach

In this section, we develop the proposed new optimization method called MAS.
Our goal is to propose a strategy that automatically combines the advantages
of an adaptive method like ADAM, with the advantages of SGD, throughout the
entire learning process. This combination of optimizers is summed as shown in
Fig. 2 where $w_{1}$ and $w_{2}$ represent the displacements on the ADAM and
SGD on the surface of the loss function, while $w1+w2$ represents the
displacement obtained thanks to our optimizer. Below we explain each line of
the MAS algorithm represented in Alg. 3.

The MAS optimizer has only two hyper-parameters which are $\lambda_{a}$ and
$\lambda_{s}$ used to balance the contribution of ADAM and SGD respectively.
In our experiments, we use only one learning rate $\eta$ for both ADAM and
SGD, but it is still possible to differentiate between the two learning rates.
In addition to the hyper-parameters typical of the MAS optimizer, all the
hyper-parameters of SGD and ADAM are also needed. In this paper, we assume to
use the most common implementation of gradient descent used in the field of
deep learning, namely the mini-batch gradient descent which divides the
training dataset into small batches that are used to calculate the model error
and update the model coefficients $w_{k}$. For each mini-batch, we calculate
the contribution derived from the two components ADAM and SGD and then update
all the coefficients as described in the three following subsections.

###  4.1 ADAM component

The complete ADAM algorithm is defined in Alg. 2. In order to use ADAM in our
optimizer, we have extracted the $\Delta_{ADAM}$ function which calculates and
returns the increments $d_{k}$ for the coefficients $w_{k}$, as defined in Eq.
16.

$d_{k}=\frac{\sqrt{1-\beta_{2}}\cdot m_{k}}{\sqrt{\widehat{v}_{k}}+\epsilon}$
(16)

The same $\Delta_{ADAM}$ function also returns the new learning rate
$\eta_{a}$ defined in Eq. 17, useful when a variable learning rate is used. In
this last case, MAS uses $\eta_{a}$ to calculate a new learning rate at each
step.

$\eta_{a}=\frac{\eta}{1-\beta_{1}}$ (17)

Now, having $\eta_{a}$ and $d_{k}$, we can directly modify the weights $w_{k}$
exactly as done in the ADAM optimizer and described in Eq. 18.

$w_{k+1}=w_{k}-\eta_{a}\cdot d_{k}$ (18)

However, we just skip this last step and use $eta_{a}$ and $d_{k}$ for our MAS
optimizer.

###  4.2 SGD component

As for the ADAM component, also the SGD component, defined in Alg. 1, has been
divided into two parts: the $\Delta_{SGD}$ function which returns the
increment to be given to the weight $w_{k}$, and the formula to update the
weights as defined in Eq. 10. The $v_{n_{k}}$ value returned by the
$\Delta_{SGD}$ function is exactly the value defined in Eq. 9, which we will
use directly for our MAS optimizer.

###  4.3 The MAS optimizer

The proposed approach can be summarized with the following Eq. 19

$w_{k+1}=w_{k}-(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot(\lambda_{s}\cdot
v_{n_{k}}+\lambda_{a}\cdot d_{k})$ (19)

where $\lambda_{s}$ is a scalar for the SGD component and $\lambda_{a}$ is
another scalar for the ADAM component used for balancing the two contribution
of the two optimizers. $\eta$ is the learning rate of the proposed MAS
optimizer, while $\eta_{a}$ is the learning rate of ADAM defined in Eq. 17.
$d_{k}$ and $v_{n_{k}}$ are the two increments define in Eq. 16 and Eq. 9
respectively.

Eq. 19 can be expanded in the following Eq. 20 to make explicit what are the
elements involved in the weights update step used by our MAS optimizer.

$\begin{split}w_{k+1}=w_{k}-(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\frac{\eta}{1-\beta_{1}})\cdot\\\
\cdot(\lambda_{s}\cdot v_{n_{k}}+\lambda_{a}\cdot\frac{\sqrt{1-\beta_{2}}\cdot
m_{k}}{\sqrt{\widehat{v}_{k}}+\epsilon})\end{split}$ (20)

where $\beta_{1}$ and $\beta_{2}$ are two parameters of the ADAM optimizer,
$v_{k}^{a}$ is defined in Eq. 12, and $m_{k}$ is defined in Eq. 11.

Algorithm 3 Mixing ADAM and SGD (MAS)

Input: the weights $w_{k}$, $\lambda_{a}$, $\lambda_{s}$, learing rate $\eta$,
weight decay $\gamma$, other SGD and ADAM paramiters …

1:for batches do

2: $d_{k},\eta_{a}=\Delta_{\text{ADAM}}(w_{k},\nabla,\eta,\gamma,\dots)$

3: $v_{n_{k}}=\Delta_{\text{SGD}}(w_{k},\nabla,\gamma,\dots)$

4: $merged=\lambda_{s}\cdot v_{n_{k}}+\lambda_{a}\cdot d_{k}$

5: $\eta_{m}=\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a}$

6: $w_{k+1}=w_{k}-\eta_{m}\cdot merged$

7:end for

The MAS algorithm can be easily implemented by following the pseudo code
defined in Alg. 3 and by calling the two functions $\Delta_{ADAM}$ defined in
Alg. 2 and $\Delta_{SGD}$ defined in Alg. 1. We can also show that convergence
is guaranteed for the MAS optimizer if we assume that convergence has been
guaranteed for the two optimizers SGD and ADAM.

######  Theorem 4.1 (MAS Cauchy necessary convergence condition).

If ADAM and SGD are two optimizers whose convergence is guaranteed then the
Cauchy necessary convergence condition is true also for MAS.

###### Proof.

Under the conditions in which the convergence of ADAM and SGD is guaranteed
[23, 24], we can say that $\sum_{k=0}^{p}{\eta\cdot v_{n_{k}}}$ and
$\sum_{k=0}^{p}{\eta_{a}\cdot d_{k}}$ converge at $\infty$ . That imply the
following:

$\lim\limits_{p\to\infty}\eta\cdot
v_{n_{p}}=\lim\limits_{p\to\infty}\eta_{a}\cdot d_{p}=0$ (21)

We can observe that
$\lim\limits_{p\to\infty}\sum_{k=0}^{p}\eta=\lim\limits_{p\to\infty}\sum_{k=0}^{p}\eta_{a}=\infty$
so we can obtain the following:

$\lim\limits_{p\to\infty}v_{n_{p}}=\lim\limits_{p\to\infty}d_{p}=0$ (22)

The thesis is that
$\sum_{k=0}^{p}(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot(\lambda_{s}\cdot
v_{n_{k}}+\lambda_{a}\cdot d_{k})$ respect the Cauchy necessary convergence
condition, so:
$\lim\limits_{p\to\infty}(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot(\lambda_{s}\cdot
v_{n_{p}}+\lambda_{a}\cdot d_{p})=0$ and for Eq. 22, this last equality is
trivially true:

$\displaystyle\lim\limits_{p\to\infty}(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot(\lambda_{s}\cdot
v_{n_{p}}+\lambda_{a}\cdot d_{p})=$
$\displaystyle(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot\lim\limits_{p\to\infty}(\lambda_{s}\cdot
0+\lambda_{a}\cdot 0)=0$ (23)

∎

######  Theorem 4.2.

(MAS convercence) If for $p\to\infty$ is valid that $\sum_{k=0}^{p}{\eta\cdot
v_{n_{k}}}=\eta\cdot m_{1}$ and $\sum_{k=0}^{p}{\eta_{a}\cdot
d_{k}}=\eta_{a}\cdot m_{2}$ where $m_{1}\in\mathbb{R}$ and
$m_{2}\in\mathbb{R}$ are two finite real values, then
$MAS=\sum_{k=0}^{p}(\lambda_{s}\cdot\eta+\lambda_{a}\cdot\eta_{a})\cdot(\lambda_{s}\cdot
v_{n_{k}}+\lambda_{a}\cdot d_{k})=\lambda_{s}^{2}\cdot\eta\cdot
m_{1}+\lambda_{a}^{2}\cdot\eta_{a}\cdot
m_{2}+\lambda_{s}\cdot\lambda_{a}\cdot\eta_{a}\cdot
m_{1}+\lambda_{s}\cdot\lambda_{a}\cdot\eta\cdot m_{2}$

###### Proof.

We can write MAS series as:

$\begin{split}MAS=\lambda_{s}^{2}\cdot\eta\cdot\sum_{k=0}^{p}v_{n_{k}}+\lambda_{a}^{2}\cdot\eta_{a}\cdot\sum_{k=0}^{p}d_{k}+\\\
+\lambda_{s}\cdot\lambda_{a}\cdot\eta_{a}\cdot\sum_{k=0}^{p}v_{n_{k}}+\lambda_{s}\cdot\lambda_{a}\cdot\eta\cdot\sum_{k=0}^{p}d_{k}\end{split}$
(24)

This can be rewritten for $p\to\infty$ as:

$MAS=\lambda_{s}^{2}\cdot\eta\cdot m_{1}+\lambda_{a}^{2}\cdot\eta_{a}\cdot
m_{2}+\lambda_{s}\cdot\lambda_{a}\cdot\eta_{a}\cdot
m_{1}+\lambda_{s}\cdot\lambda_{a}\cdot\eta\cdot m_{2}$ (25)

∎

Figure 3:  Behavior of the three optimizers MAS, ADAM and SGD on the surface
defined in Eq.28. For better visualization the SGD was shifted on X axis of
$0.1$.

###  4.4 Geometric explanation

We can see optimizers as two explorers $w_{1}$ and $w_{2}$ who want to explore
an environment (the surface of a loss function). If the two explorers agree to
go in a similar direction, then they quickly go in that direction
($w_{1}+w_{2}$). Otherwise, if they disagree and each prefers a different
direction than the other, then they proceed more cautiously and slower
($w_{1}+w_{2}$). As we can see in Fig. 2a, if the directions of the
displacement of $w_{1}$ and $w_{2}$ are similar then the amplitude of the
resulting new displacement $w_{1}+w_{2}$ is increased, however, as shown in
Fig. 2b, if the directions of the two displacements $w_{1}$ and $w_{2}$ are
not similar then the amplitude of the new displacement $w_{1}+w_{2}$ has
decreased.

In our approach, the sum $w_{1}+w_{2}$ is weighted (see red vectors in Fig.
2a) so one of the two optimizers SGD or ADAM can become more relevant than the
other in the choice of direction for MAS, hence the direction resultant may
tend towards one of the two. In MAS we set the weight of the two contributions
so as to have a sum $\lambda_{1}+\lambda_{2}=1$ in order to maintain a
learning rate of the same order of magnitude.

Another important component that greatly affects the MAS shift module at each
training step is its learning rate defined in Eq. 19 which combines $\eta$ and
$\eta_{a}$. The shifts are scaled using the learning rate, so there is a
situation where MAS gets more thrust than the ADAM and SGD starting shifts. In
particular, we can imagine that the displacement vector of ADAM has a greater
magnitude than SGD and the learning rate of SGD is greater than that of ADAM.
In this case, the MAS shift has a greater vector magnitude than SGD and a
higher ADAM learning rate which can cause a large increase in the MAS shift
towards the search of a minimum.

Figure 4:  Behavior of the three optimizers MAS, ADAM and SGD on the
Rosenbrook’s surface with $a=1$ and $b=100$ Figure 5:  Behavior of the three
optimizers MAS, ADAM and SGD on the surface $z=\frac{\lvert
x\lvert}{10}+\lvert y\lvert$

###  4.5 Toy examples

To better understand our proposal, we built a toy example where we highlight
the main behaviour of MAS. More precisely we consider the following example:

$x=[1,2],~{}y=[2,4]$ (26) $p_{i}=w_{1}\cdot(w_{2}\cdot x_{i})$ (27)
$\mathcal{L}(w_{1},w_{2})=\sum_{i=0}^{1}(p_{i}-y_{i}^{2})$ (28)

We set $\beta_{1}=0.9$, $\beta_{2}=0.999$, $\epsilon=10^{-8}$,
$amsgrad=False$, dampening $d=0$, $nesterov=False$ and $\mu=0$. As we can see
in Fig. 3 our MAS optimizer goes faster towards the minimum value after only
two epochs, SGD is fast at the first epoch, however, it decreases its speed
soon after and comes close to the minimum after 100 epochs, ADAM instead
reaches its minimum after 25 epochs. Our approach can be fast when it gets a
large $v_{k}$ from SGD and a large $\eta_{a}$ from ADAM.

Another toy example can be done with the benchmark Rosenbrook [25] function:

$z=(a-y)^{2}+b\cdot(y-x^{2})^{2}$ (29)

We set $a=1$ and $b=100$, weight $x=3$ and weight $y=1$, $lr=0.0001$,
$epochs=1000$, and default paramiter for ADAM and SGD. The MAS optimizer sets
$\lambda_{s}=\lambda_{a}=0.5$. The comparative result for the minimization of
this function is shown in Fig. 4. In this experiment, we can see how by
combining the two optimizers ADAM and SGD we can obtain a better result than
the single optimizers. For this function, going from the starting point
towards the direction of the maximum slope means moving away from the minimum,
and therefore it takes a lot of training epochs to approach the minimum.

Let’s use a final toy example to highlight the behavior of the MAS optimizer.
In this case we look for the minimum of the function $z=\frac{\lvert
x\lvert}{10}+\lvert y\lvert$. We set the weights $x=3$ and $y=2$, $lr=0.01$,
$epochs=400$ and use all the default parameters for ADAM and SGD. MAS assigns
the same value $0.5$ for the two lamdas hyper-parameters. In Fig. 5 we can see
how not all the paths between the paths of ADAM and SGD are the best choice.
Since MAS, as shown in Fig.2, goes towards an average direction with respect
to that of ADAM and SGD, then in this case ADAM arrives first at the minimum
point.

Figure 6: On the left some sample images for each of the 10 classes of
CIFAR-10 (one class for each row). On the right 10 classes randomly selected
from the set of 100 classes of CIFAR-100.

##  5 Datasets

In this section, we briefly describe the datasets used in the experimental
phase.

The Cifar10 [26] dataset consists of 60,000 images divided into 10 classes
(6000 per class) with a training set size and test set size of 50000 and 10000
respectively. Each input sample is a low-resolution color image of size
$32\times 32$. The 10 classes are airplanes, cars, birds, cats, deer, dogs,
frogs, horses, ships, and trucks.

Table 1: Accuracy results on Cifar10, after 6 runs of 200 epochs. Name | $\lambda_{a}$ | $\lambda_{s}$ | avg. acc | acc max
---|---|---|---|---
Resnet18 |  |  |  |
Adam | 1 | 0 | 84.68 | 86.24
SGD | 0 | 1 | 78.87 | 79.19
MAS | 0.5 | 0.5 | 85.36 | 85.80
MAS | 0.4 | 0.6 | 85.64 | 86.59
MAS | 0.6 | 0.4 | 85.89 | 86.56
MAS | 0.7 | 0.3 | 85.39 | 86.09
MAS | 0.3 | 0.7 | 85.57 | 86.85
Resnet34 |  |  |  |
Adam | 1 | 0 | 82.98 | 83.52
SGD | 0 | 1 | 82.92 | 83.25
MAS | 0.5 | 0.5 | 84.99 | 85.69
MAS | 0.4 | 0.6 | 85.75 | 86.12
MAS | 0.6 | 0.4 | 84.63 | 85.27
MAS | 0.7 | 0.3 | 84.46 | 84.80
MAS | 0.3 | 0.7 | 85.71 | 86.14

The Cifar100 [26] dataset consist of 60000 images divided in 100 classes (600
per classes) with a training set size and test set size of 50000 and 10000
respectively. Each input sample is a $32\times 32$ colour images with a low
resolution. In Fig. 6 we report some representative examples of the two
datasets, extracted from a subset of classes.

The Corpus of Linguistic Acceptability (CoLA) [27] is another dataset which
contains 9594 sentences belonging to training and validation sets, and
excludes 1063 sentences belonging to a set of tests kept out. In our
experiment, we only used the training set and the test set.

The AG’s news corpus [28, 29] is the last dataset used in our experiments. It
is a dataset that contains news articles from the web subdivided into four
classes. It has 30,000 training samples and 1900 test samples.

Table 2: Accuracy results on Cifar100, after 7 runs of 200 epochs. Name | $\lambda_{a}$ | $\lambda_{s}$ | avg acc. | acc max
---|---|---|---|---
Resnet18 |  |  |  |
Adam | 1 | 0 | 49.56 | 50.28
SGD | 0 | 1 | 49.48 | 50.43
MAS | 0.5 | 0.5 | 55.08 | 56.68
MAS | 0.4 | 0.6 | 56.23 | 56.83
MAS | 0.6 | 0.4 | 53.82 | 54.44
MAS | 0.7 | 0.3 | 52.92 | 54.07
MAS | 0.3 | 0.7 | 58.01 | 58.48
Resnet34 |  |  |  |
Adam | 1 | 0 | 50.66 | 51.92
SGD | 0 | 1 | 52.81 | 53.45
MAS | 0.5 | 0.5 | 51.52 | 53.26
MAS | 0.4 | 0.6 | 51.73 | 53.48
MAS | 0.6 | 0.4 | 52.15 | 53.95
MAS | 0.7 | 0.3 | 53.06 | 54.50
MAS | 0.3 | 0.7 | 51.96 | 53.32

##  6 Experiments

The optimizer MAS proposed is a generic solution not oriented exclusively to
image analysis, so we conduct experiments on both image classification and
text document classification. By doing so, we are able to give a clear
indication of the behavior of the proposed optimizer in different contexts,
also bearing in mind that many problems, such as audio recognition, can be
traced back to image analysis. In all the experiments $\beta_{1}=0.9$,
$\beta_{2}=0.999$, $\epsilon=10^{-8}$, $amsgrad=False$, dampening $d=0$ and
$nesterov=False$.

###  6.1 Experiments with images

In this first group of experiments we use two well-known image datasets for:
(1) conduct an analysis of the two main parameters of MAS, $\lambda_{a}$ and
$\lambda_{s}$; (2) compare the performance of MAS with respect to the two
starting optimizers SGD and ADAM; (3) analyze the behavior of MAS with
different neural models.

The datasets used in this first group of experiments are Cifar10 and Cifar100.
The neural models compared are two ResNet and in particular, we use Resnet18
and Resnet34 [30, 31]. We analyzed $\lambda_{a}$ and $\lambda_{s}$ when they
assume values from the set $\\{0.3,0.4,0.5,0.6,0.7\\}$ so that
$\lambda_{a}+\lambda_{s}=1$. The numerical results are all grouped in the two
Tabs. 1 and 2.

For Cifar10 we report the average accuracies calculated on 6 runs of each
experiment and for 200 epochs. We set $\eta=0.001$, momentum $\mu=0.95$, and
batch size equal to 1024. From the results reported in Tab. 1, we can observe
that the best average results are obtained in correspondence with
$\lambda_{s}=0.6$ and $\lambda_{a}=0.4$ but we can also note that all the
other settings used for MAS allow obtaining better results than ADAM and SGD
for both Resnet18 and Resnet34. To better understand what happens during the
training phase, in Fig. 7 we represent the accuracy of the test and the
corresponding loss values of the experiments that produced the best results
with Resnet18 on Cifar10. As we can see even though the best execution of ADAM
has the lowest loss value, our optimizer with $\lambda_{a}=0.3$ and
$\lambda_{s}=0.7$ offers better testing accuracy. The other combination
represented with $\lambda_{a}=0.5$ and $\lambda_{s}=0.5$ is also very similar
to the better accuracy of ADAM and in any case always better than the average
accuracy produced by ADAM. So in general we can say that the MAS optimizer
leads to a better generalize than the other optimizers used.

Figure 7: Resnet18 test accuracies and test loss of the best results obtained
on Cifar10.

For Cifar100 we use the same settings as the Cifar10 experiment except for
$\eta=0.008$ and batch size set to 512. We can see the results in Tab. 2. Also
for this dataset, we can see that for Resnet18 all MAS configurations work
better than SGD and ADAM. Instead, looking at the results obtained with the
Resnet34, we can say that only one configuration of MAS exceeds the average
accuracies of SGD and ADAM, but if we look at the maximum accuracy values,
more than one configuration of MAS is better than the results available with
ADAM and SGD.

In conclusion, as we have seen from the results shown in Tab. 1 and Tab. 2,
the proposed method leads to a better generalization than the other optimizers
used in each experiment. We get better results both by setting $\lambda_{a}$
and $\lambda_{s}$ well, and also even when we don’t use the best set of
parameters.

###  6.2 Experiments with text documents

In this last group of experiments, we use the two datasets of text documents:
CoLA and AG’s News. As neural model, we use a model based on BERT [17] which
is one of the best techniques for working with text documents. To run fewer
epochs we use a pre-trained version [32] of BERT. In these experiments, we
also introduce the comparison with the AdamW optimizer which is usually the
optimizer used in BERT-based models.

For the CoLA dataset we set $\eta=0.0002$, momentum $\mu=0.95$, and batch size
equal to 100. We ran the experiments 5 times for 50 epochs. For the AG’s News
dataset we set the same parameters used for CoLA, but we only run it for 10
epochs because it gets good results in the firsts epochs and also because the
dataset is very large and therefore takes more time.

We can see all the results in Tab. 3. Even for text analysis problems, we can
confirm the results of the experiments done on images: although AdamW
sometimes has better performances than ADAM, our proposed optimizer performs
better than other optimizers used in this paper.

Table 3: Accuracy results of BERT pre-trained on CoLA (50 epochs) and AG’s news (10 epochs), after 5 runs. Name | $\lambda_{a}$ | $\lambda_{s}$ | avg acc. | acc max
---|---|---|---|---
CoLA |  |  |  |
AdamW | - | - | 78.59 | 85.96
Adam | 1 | 0 | 79.85 | 83.30
SGD | 0 | 1 | 81.48 | 81.78
MAS | 0.5 | 0.5 | 85.92 | 86.72
MAS | 0.4 | 0.6 | 86.18 | 87.66
MAS | 0.6 | 0.4 | 85.45 | 86.34
MAS | 0.7 | 0.3 | 84.66 | 85.78
MAS | 0.3 | 0.7 | 86.34 | 86.91
AG’s News |  |  |  |
AdamW | - | - | 92.62 | 92.93
Adam | 1 | 0 | 92.55 | 92.67
SGD | 0 | 1 | 91.28 | 91.39
MAS | 0.5 | 0.5 | 93.72 | 93.80
MAS | 0.4 | 0.6 | 93.82 | 93.98
MAS | 0.6 | 0.4 | 93.55 | 93.67
MAS | 0.7 | 0.3 | 93.19 | 93.32
MAS | 0.3 | 0.7 | 93.86 | 93.99

##  7 Conclusion

In this paper, we introduced MAS (Mixing ADAM and SGD) a new Combined
Optimization Method that combines the capability of two different optimizers
into one. We demonstrate by experiments the capability of our proposal to
overcome the single optimizers used in our experiments and achieve better
performance. To balance the contribution of the optimizers used within MAS, we
introduce two new hyperparameters $\lambda_{a}$, $\lambda_{s}$ and show
experimentally that in almost all configurations of these parameters, the
results are better than the results obtained with the other single optimizers.
In future work, it is possible to change ADAM and SGD and try to mix different
optimizers also without the limitations of using only two optimizers. Another
significant future work is to try to change dynamically during the training
the influence (lambda hyper-parameters) of the two combined optimizers, in the
hypothesis that this can improve further the generalization performance.

## References

  * [1] H. Robbins and S. Monro, “A stochastic approximation method,” The annals of mathematical statistics, pp. 400–407, 1951.
  * [2] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [3] M. Zaheer, S. Reddi, D. Sachan, S. Kale, and S. Kumar, “Adaptive methods for nonconvex optimization,” in Advances in neural information processing systems, pp. 9793–9803, 2018.
  * [4] L. Luo, Y. Xiong, Y. Liu, and X. Sun, “Adaptive gradient methods with dynamic bound of learning rate,” arXiv preprint arXiv:1902.09843, 2019.
  * [5] S. Bera and V. K. Shrivastava, “Analysis of various optimizers on deep convolutional neural network model in the application of hyperspectral remote sensing image classification,” International Journal of Remote Sensing, vol. 41, no. 7, pp. 2664–2683, 2020.
  * [6] A. Graves, “Generating sequences with recurrent neural networks,” arXiv preprint arXiv:1308.0850, 2013.
  * [7] J. Duchi, E. Hazan, and Y. Singer, “Adaptive subgradient methods for online learning and stochastic optimization.,” Journal of machine learning research, vol. 12, no. 7, 2011.
  * [8] M. D. Zeiler, “Adadelta: an adaptive learning rate method,” arXiv preprint arXiv:1212.5701, 2012.
  * [9] L. Luo, Y. Xiong, Y. Liu, and X. Sun, “Adaptive gradient methods with dynamic bound of learning rate,” in International Conference on Learning Representations, 2018.
  * [10] T. Kobayashi, “Scw-sgd: Stochastically confidence-weighted sgd,” in 2020 IEEE International Conference on Image Processing (ICIP), pp. 1746–1750, IEEE, 2020.
  * [11] Z. Zhang, “Improved adam optimizer for deep neural networks,” in 2018 IEEE/ACM 26th International Symposium on Quality of Service (IWQoS), pp. 1–2, IEEE, 2018.
  * [12] N. Landro, “Mixing adam and sgd: a combined optimization method with pythorch,” 2020.  https://gitlab.com/nicolalandro/multi_optimizer.
  * [13] N. S. Keskar and R. Socher, “Improving generalization performance by switching from adam to sgd,” arXiv preprint arXiv:1712.07628, 2017.
  * [14] X. Cui, W. Zhang, Z. Tüske, and M. Picheny, “Evolutionary stochastic gradient descent for optimization of deep neural networks,” in Advances in neural information processing systems, pp. 6048–6058, 2018.
  * [15] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” arXiv preprint arXiv:1711.05101, 2017.
  * [16] I. Loshchilov and F. Hutter, “Fixing weight decay regularization in adam,”
, 2018.

  * [17] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” arXiv preprint arXiv:1810.04805, 2018.
  * [18] A. Krogh and J. A. Hertz, “A simple weight decay can improve generalization,” in Advances in neural information processing systems, pp. 950–957, 1992\.
  * [19] I. Sutskever, J. Martens, G. Dahl, and G. Hinton, “On the importance of initialization and momentum in deep learning,” in International conference on machine learning, pp. 1139–1147, 2013.
  * [20] G. Damaskinos, E. M. E. Mhamdi, R. Guerraoui, R. Patra, and M. Taziki, “Asynchronous byzantine machine learning (the case of sgd),” arXiv preprint arXiv:1802.07928, 2018.
  * [21] C. Liu and M. Belkin, “Accelerating sgd with momentum for over-parameterized learning,” arXiv preprint arXiv:1810.13395, 2018.
  * [22] J. Chen, D. Zhou, Y. Tang, Z. Yang, and Q. Gu, “Closing the generalization gap of adaptive gradient methods in training deep neural networks,” arXiv preprint arXiv:1806.06763, 2018.
  * [23] S. J. Reddi, S. Kale, and S. Kumar, “On the convergence of adam and beyond,” arXiv preprint arXiv:1904.09237, 2019.
  * [24] J. D. Lee, M. Simchowitz, M. I. Jordan, and B. Recht, “Gradient descent only converges to minimizers,” in Conference on learning theory, pp. 1246–1257, 2016.
  * [25] H. Rosenbrock, “An automatic method for finding the greatest or least value of a function,” The Computer Journal, vol. 3, no. 3, pp. 175–184, 1960.
  * [26] A. Krizhevsky, G. Hinton, et al., “Learning multiple layers of features from tiny images,” Citeseer, 2009.
  * [27] A. Warstadt, A. Singh, and S. R. Bowman, “Neural network acceptability judgments,” arXiv preprint arXiv:1805.12471, 2018.
  * [28] A. Gulli, “Ag’s corpus of news articles,” 2005.  http://groups.di.unipi.it/ gulli/
AG_corpus_of_news_articles.html.

  * [29] X. Zhang, J. Zhao, and Y. LeCun, “Character-level convolutional networks for text classification,” in Advances in neural information processing systems, pp. 649–657, 2015.
  * [30] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” arXiv preprint arXiv:1512.03385, 2015.
  * [31] S. Targ, D. Almeida, and K. Lyman, “Resnet in resnet: Generalizing residual architectures,” arXiv preprint arXiv:1603.08029, 2016.
  * [32] huggingface.co, “Bert base uncased pre-trained model.”  https://huggingface.co/bert-base-uncased.