# Nonmagnetic Ground State in RuO2 Revealed by Muon Spin Rotation

M. Hiraishi<EMAIL_ADDRESS>Graduate School of
Science and Engineering, Ibaraki University, Mito, Ibaraki 310-8512, Japan
Muon Science Laboratory, Institute of Materials Structure Science, High Energy
Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan  H.
Okabe  Institute for Materials Research, Tohoku University (IMR), Aoba-ku,
Sendai 980-8577, Japan  Muon Science Laboratory, Institute of Materials
Structure Science, High Energy Accelerator Research Organization (KEK),
Tsukuba, Ibaraki 305-0801, Japan  A. Koda  Muon Science Laboratory, Institute
of Materials Structure Science, High Energy Accelerator Research Organization
(KEK), Tsukuba, Ibaraki 305-0801, Japan  Graduate University for Advanced
Studies, SOKENDAI  R. Kadono  Muon Science Laboratory, Institute of Materials
Structure Science, High Energy Accelerator Research Organization (KEK),
Tsukuba, Ibaraki 305-0801, Japan  T. Muroi  Institute for Solid State Physics,
University of Tokyo, Kashiwa, Chiba 277-8581, Japan  D. Hirai  Department of
Applied Physics, Graduate School of Engineering, Nagoya University, Chikusa-
ku, Nagoya 464-8603, Japan  Z. Hiroi  Institute for Solid State Physics,
University of Tokyo, Kashiwa, Chiba 277-8581, Japan

###### Abstract

The magnetic ground state of single crystalline RuO2 was investigated by the
muon spin rotation/relaxation ($\mu$SR) experiment. The spin precession signal
due to the spontaneous internal magnetic field $B_{\rm loc}$, which is
expected in the magnetically ordered phase, was not observed in the
temperature range 5–400 K. Muon sites were evaluated by first-principles
calculations using dilute hydrogen simulating muon as pseudo-hydrogen, and
$B_{\rm loc}$ was simulated for the antiferromagnetic (AFM) structures with a
Ru magnetic moment $|{\bm{m}}_{\rm Ru}|\approx 0.05\mu_{\mathrm{B}}$ suggested
from diffraction experiments. As a result, the possibility was ruled out that
muons are localized at sites where $B_{\rm loc}$ accidentally cancels.
Conversely, assuming that the slow relaxation observed in $\mu$SR spectra was
part of the precession signal, the upper limit for the magnitude of
$|{\bm{m}}_{\rm Ru}|$ was estimated to be $4.8(2)\times 10^{-4}\mu_{\rm B}$,
which is significantly less than $0.05\mu_{\mathrm{B}}$. These results
indicate that the AFM order, as reported, is unlikely to exist in the bulk
crystal.

Ruthenium oxide RuO2 with rutile structure is a material that has been applied
in a variety of fields due to its high catalytic activity and remarkable
chemical stability Over (2012). From the viewpoint of electronic properties,
RuO2 has been considered to be an ordinary Pauli paramagnetic (i.e.,
nonmagnetic) metal Ryden and Lawson (2003). However, it has recently been
reported to have a topological electronic structure Sun  _et al._ (2017);
Jovic  _et al._ (2018). Furthermore, based on neutron diffraction Berlijn  _et
al._ (2017) and resonant X-ray scattering experiments Zhu  _et al._ (2019),
RuO2 has been claimed to exhibit antiferromagnetic (AFM) ordering with a high
Néel temperature ($>300$ K) and a Ru magnetic moment size of
$\sim$$0.05\mu_{\mathrm{B}}$. The presence of AFM phase has generated interest
in its potential application for spintronic devices Wadley  _et al._ (2016);
Baltz  _et al._ (2018). In addition, strain-induced superconductivity has also
been reported Uchida  _et al._ (2020); Ruf  _et al._ (2021), thereby fueling
interest in its detailed electronic properties.

Following the inference of AFM ordering in RuO2, theoretical predictions and
experimental results for various anomalies related to transport phenomena have
been reported. Anomalous Hall effects associated with collinear AFM phase and
the non-concentric symmetric position of the nonmagnetic oxygen have been
predicted from theoretical studies $\mathrm{\check{S}}$mejkal  _et al._
(2020), and experimental result supporting this prediction has been reported
Feng  _et al._ (2022). In addition, spin current due to the spin-splitter
effect generated in the AFM phase has been theoretically proposed González-
Hernández  _et al._ (2021), which has been followed by reports of experimental
results in favor of the prediction Bai  _et al._ (2022); Bose  _et al._
(2022); Karube  _et al._ (2022). Occurrence of the chirality magneto-optical
effect is also theoretically predicted Zhou  _et al._ (2021).

However, the reported size of Ru magnetic moments is close to the limit of
sensitivity in neutron and X-ray diffraction experiments. Moreover, a recent
theoretical study suggests that AFM ordering may be induced by hole doping due
to Ru vacancies in RuO2 which is intrinsically nonmagnetic Smolyanyuk  _et
al._ (2023). Thus, verification of the AFM phase with local magnetic probes
that are complementary to diffraction experiment is highly required. Motivated
by this situation, we have investigated the magnetic ground state of RuO2 by
the muon spin rotation/relaxation ($\mu$SR). The experimental results show
that the spontaneous muon spin precession signal expected in the magnetically
ordered phase is not observed in the temperature range of 5–400 K. The first-
principles density functional theory (DFT) calculations of muon sites exclude
the possibility that muons are localized at sites where the internal magnetic
field exerted by the Ru atoms is canceled by chance. These results support the
scenario that the bulk crystal RuO2 is a nonmagnetic metal.

$\mu$SR is a magnetically sensitive probe in which spin-polarized muons
($\mu^{+}$) stopped in the sample directly infer the internal magnetic field
$B_{\rm loc}$ at the interstitial sites via the spontaneous Larmor precession.
Muons are provided as a 100%-spin-polarized beam by proton accelerator
facilities, which enables $\mu$SR measurements in zero magnetic field. The
sensitivity of $\mu$SR as a local probe has been demonstrated in the discovery
of AFM ordering in the parent compounds of high-$T_{\mathrm{c}}$ cuprate
superconductors Nishida  _et al._ (1987); Uemura  _et al._ (1987); Luke  _et
al._ (1989), as well as in the study of superconductivity Aeppli  _et al._
(1987); Uemura  _et al._ (1989); Sonier  _et al._ (2000). In addition, $\mu$SR
has various advantages: the muon implantation energy is high enough ($\geq 4$
MeV) to be surface independent (bulk sensitive), and the implanted muons
(volume concentration $\sim$$10^{5}$ cm-3) decay with an average lifetime of
2.2 µs, so they do not accumulate in the sample unlike other ion beam
irradiation measurements. Therefore, the detection of $B_{\rm loc}$ due to
ordered Ru moments by $\mu$SR could provide direct experimental evidence for
the occurrence of AFM ordering.

For the quantitative assessment of $B_{\rm loc}$ using $\mu$SR, it is
necessary to obtain information on the local electronic structure of $\mu^{+}$
in RuO2; it behaves as a light isotope of proton/hydrogen (pseudo-H, hereafter
denoted by the element symbol “Mu”) in matter. To this end, first-principles
DFT calculations were performed using Quantum Espresso Giannozzi  _et al._
(2009, 2017, 2020) to investigate the Mu stopping sites (DFT+Mu). Structural
relaxation calculations with H (to mimic Mu) using the GGA-PBE exchange
correlation function were performed on a $3\times 3\times 3$ superlattice (54
Ru, 108 oxygen and H) of a rutile structure ($P4_{2}/mnm$) as the initial
structure with cutoff energy of 60 Ry and $K$ points were set to $3\times
3\times 4$. Structure was relaxed until the maximum force on each atom was
less than $1\times 10^{-3}$ (Ry/Bohr) in the nonmagnetic state. For the on-
site Coulomb interaction $U$ between Ru $4d$ orbitals, calculations were
performed for values of $U=0$–3 eV and results were found to be mostly
independent of $U$. We then used the result for $U=2.0$ eV (where the AFM
state is theoretically favored Berlijn  _et al._ (2017);
$\mathrm{\check{S}}$mejkal  _et al._ (2020); Ruf  _et al._ (2021)) for the
simulation of $B_{\rm loc}$. More detailed information can be found in the
Supplemental Material (SM) sm , which includes Refs. RodrÃguez-Carvajal
(1993); Dolomanov  _et al._ (2009); Bourhis  _et al._ (2015); Sheldrick
(2015); PP ; Hamann (2013).

Single crystals of RuO2 were grown from a polycrystalline sample in oxygen
flow. RuO2 reagent (Rare Metallic Co., Ltd., 99.95%) was placed in the high
temperature section of the furnace (1100∘C) and held for two weeks while
oxygen flowed at a rate of 30 ml/min to obtain blue-black metallic single
crystals in the low-temperature section downstream of the furnace. The
crystals were about $3\times 3\times 0.5$ mm3 at the maximum with a thin
rhombic columnar shape. To evaluate the quality of these crystals, single-
crystal (sc) XRD analysis and electrical resistivity measurements were
performed on the samples from the same batch (for more details, see SM sm ).
The Ru occupancy was deduced to be 1.02(2) from sc XRD, indicating that the
crystal is free of Ru defects at the percent level. The resistivity
measurements showed a maximum residual resistivity ratio (rrr) of $\sim$1500\.
This is much higher than previously reported values ($\sim$30 Ruf  _et al._
(2021) or $\sim$200 Rogers  _et al._ (1969)), confirming that the crystals
were of extremely high quality. These crystals were loaded on a silver holder
covering about $15\times 15$ mm2 in a mosaic pattern (corresponding to the
powder average for the $\mu$SR spectra). As a reference, $\mu$SR measurements
on RuO2 powder (same supplier, 99.9%) as delivered were performed, and
obtained qualitatively similar results to those of the single-crystalline
samples. However, the results showed certain differences in details presumably
attributed to the sample qualities such as impurities and defects, whose
details are found in SM sm .

Conventional $\mu$SR measurements were performed in zero field (ZF), weak
longitudinal field (LF) and weak transverse field (TF) using the S1 instrument
(ARTEMIS) Kojima  _et al._ (2014), a general-purpose $\mu$SR spectrometer
installed at the S1 area of the Materials and Life science Experiment Facility
(MLF), J-PARC. A 100% spin-polarized $\mu^{+}$ beam ($\sim$4 MeV) was stopped
in the sample mounted on the He gas-flow cryostat for controlling sample
temperature. The time-dependent $\mu^{+}$ spin polarization, which reflects
the magnetic field distribution at the Mu site, was monitored via the forward-
backward asymmetry [$A(t)$] of positrons. The $\mu$SR spectra were analyzed by
least-square curve fits using “musrfit” Suter and Wojek (2012). The background
contribution which missed the sample was calibrated using a Ho sample of the
same cross section as RuO2 to be $A_{\mathrm{BG}}=0.0815$, and subtracted from
$A(t)$. The evaluation of $B_{\rm loc}$ due to Ru magnetic moments and
putative Ru $4d$ ordered moments were performed using “dipelec” code Kojima .

Figure 1(a) shows ZF-$\mu$SR time spectra which are normalized by the
asymmetry at $t=0$ to yield time-dependent spin polarization,
$G_{z}(t)=A(t)/A_{0}$ [with $A_{0}=A(0)$]. The line shape of these spectra is
well represented by slow exponential damping over the entire temperature range
of 5–400 K, and no sinusoidal component for the Larmor precession is observed.
This indicates that there is no homogeneous $B_{\rm loc}$ established at the
Mu site. Moreover, as shown in Fig. 1(b), the slow relaxation observed in ZF
is completely suppressed by applying LF of 1 mT, indicating that
$B_{\mathrm{loc}}$ is quasistatic.

Figure 1: (a) Examples of ZF-$\mu$SR spectra at typical temperatures, and (b)
magnetic field dependence of LF-$\mu$SR spectra at 5 K observed in single-
crystalline RuO2. Solid curves represent the least-square fit by Eq. (1).
Inset: the spontaneous spin precession signal expected in the presence of the
presumed AFM order (see text).

Consequently, the curve fits of the spectra was performed using the following
function,

$\displaystyle A_{0}G_{z}(t)$
$\displaystyle=A_{1}G_{\mathrm{KT}}^{\mathrm{L}}(\lambda,B_{\mathrm{LF}},t),$
(1) $\displaystyle G_{\mathrm{KT}}^{\mathrm{L}}(\lambda,0,t)$
$\displaystyle=\frac{1}{3}+\frac{2}{3}(1-\lambda t)\exp(-\lambda t),$

where $G_{\mathrm{KT}}^{\mathrm{L}}(\lambda,B_{\mathrm{LF}},t)$ is the
relaxation function for a static-Lorentzian field distribution with $\lambda$
being the relaxation rate Uemura  _et al._ (1985), $B_{\mathrm{LF}}$ is the
external magnetic field.

The result of curve fits using Eq. (1) is shown for the ZF time spectra at
various temperatures and for the LF spectra at 5 K (including their $B_{\rm
LF}$ dependence) in Fig. 1(a) and (b), both of which show good agreement.
Considering that spin relaxation due to nuclear magnetic moments is usually
well described by Gaussian Kubo-Toyabe functions Hayano  _et al._ (1979), the
origin of the observed Lorentzian-like behavior is unknown at this stage.
Although the issue remains to be addressed in future work, these results
indicate that $\lambda$ is not subject to fluctuation by self-diffusion of Mu
nor magnetic fluctuation, ruling out the possibility that $B_{\rm loc}$ is
effectively averaged out to zero due to the motional narrowing effect (see
also the discussion below).

The temperature dependence of the initial Mu polarization [$G_{z}(0)$] is
shown in Fig. 2(a), which is nearly unity irrespective of temperature for
5–400 K, indicating that there is no muonium (a neutral H atom-like state,
Mu0) formation nor fast depolarization components beyond the time resolution
of the instrument. Thus, implanted Mu are mostly in the diamagnetic state (Mu+
or Mu-). Considering that RuO2 is a metal, Mu is assumed to be an isolated
state of Mu+ due to shielding by conduction electrons.

As shown in Fig. 2(b), $\lambda$ is nearly constant at temperatures below
about 100 K, while it decreases slowly at higher temperatures. Assuming that
this is due to thermally activated diffusive motion of Mu, the temperature
dependence of $\lambda$ is given as

$\displaystyle\lambda(T)=\lambda_{0}\left[1+N\exp(-E_{\mathrm{a}}/k_{\mathrm{B}}T)\right]^{-1},$
(2)

where $\lambda_{0}$ is the relaxation rate at $T\rightarrow 0$ K,
$E_{\mathrm{a}}$ is the activation energy, and $N$ is a physical quantities
related to the density of states Cox  _et al._ (2006). As shown by dashed line
in Fig. 2(b), curve fit provides good agreement with the data, yielding
$E_{\mathrm{a}}=62(8)$ meV being also typical for Mu diffusion in oxides Ito
_et al._ (2023). Thus, it suggests that Mu diffusive motion may occur above
$\sim$100 K, suggesting the possibility that $B_{\rm loc}$ is effectively
reduced by the motional effect. Conversely, however, it suggests that at least
below $\sim$100 K, $B_{\rm loc}$ is free from such motional effect.

Figure 2: Temperature dependence of (a) the initial polarization and (b)
relaxation rate obtained by the analysis using Eq. (1). The dashed line in (b)
is the result of a curve fit using Eq. (2).

In general, there are interstitial sites with high symmetry in
antiferromagnets where $B_{\rm loc}$ exerted from magnetic ions cancels out.
Given the reported magnetic structure Berlijn  _et al._ (2017); Zhu  _et al._
(2019), that is the case for the $8j$ site $(0.25,0.25,0.25)$ situated at the
midpoint between Ru $(0,0,0)$ and Ru $(0.5,0.5,0.5)$, and the $4d$ site
$(0,0.5,0.25)$. Provided that Mu accidentally occupies these sites, the
spontaneous spin precession signal due to AFM ordering will be absent, thus
insensitive to $B_{\rm loc}$. Therefore, it is important to gain reliable
estimation of the Mu site(s) and associated sensitivity to $B_{\rm loc}$
expected for the reported magnetic structure.

First, the most stable site (labeled Site1) for H inferred from structural
relaxation calculation is shown in Fig. 3. The H atom is about 0.1 nm from the
nearest oxygen, and the bond direction is approximately in the $ab$ plane
(there are four equivalent sites in the unit cell). This reproduces the
previously reported result Jadon  _et al._ (2021), and is also in line with
the general trend in oxides including rutile TiO2 Marinopoulos  _et al._
(2018) that Mu as pseudo-H forms OH bond with oxygen. According to the total
energy calculations performed with H fixed at the $8j$ and $4d$ sites, the
total energies are about 2.1 and 1.4 eV higher than Site1, respectively, which
rules out the possibility that Mu is stationary at these sites where $B_{\rm
loc}$ cancels out.

The mass of Mu is approximately 1/9 that of H, and quantum effects such as
zero-point energy $E_{0}$ can be pronounced. In particular, when $E_{0}$ is
greater than the potential barrier $V_{\mathrm{b}}$ of the jump diffusion path
to other stable sites, Mu may be delocalized to occupy site(s) away from
Site1: such an effect has been actually reported for the iron-based
superconductor LaFeAsO0.49H0.51 Hiraishi  _et al._ (2020). According to
earlier reports, the activation energy ($\approx V_{\mathrm{b}}$) of the
diffusion path from the Site1 to the next OMu bonding site at the same $z$ in
the same $c$-channel [$(x,y,z)\rightarrow(1-x,-y,z)$ in the unit cell] is
estimated to be 0.27 or 0.217 eV (the latter includes the zero-point energy
correction) Kim and Lai (2015); Jadon  _et al._ (2021). When $E_{0}$ of Mu
(roughly three times that of H in the harmonic approximation) is greater than
this barrier, Mu may be localized at the center in the $c$-channel. This site
is the $4c$ site $(0.5,0,0)$ in the unit cell (denoted Site2, as shown in Fig.
3).

Figure 3: Local atomic configurations of hydrogen defect in RuO2 obtained by
DFT calculations with lattice relaxation in $3\times 3\times 3$ superlattice.
Site1: the lowest energy site (forming OH bond). Site2 (light blue): the
$c$-channel center ($4d$) site. The structure was displayed using VESTA Momma
and Izumi (2011).

Considering that the direction of the Ru moment ${\bm{m}}_{\mathrm{Ru}}$ is
not uniquely provided in the literature Berlijn  _et al._ (2017); Zhu  _et
al._ (2019) and that ${\bm{m}}_{\mathrm{Ru}}$ is theoretically predicted to
vary from (001) to (100) as the electron filling progresses
$\mathrm{\check{S}}$mejkal  _et al._ (2020), $B_{\rm loc}$ is calculated for
the arbitrary direction of ${\bm{m}}_{\mathrm{Ru}}$. Our calculations reveal
that the $B_{\mathrm{loc}}$ is never canceled in any direction of
${\bm{m}}_{\mathrm{Ru}}$ (see SM sm ). In the case of
${\bm{m}}_{\mathrm{Ru}}\parallel(001)$, the corresponding $B_{\rm loc}$ is
6.17 mT for Site1 and 5.96 mT for Site2, yielding the precession frequencies
0.84 and 0.81 MHz for Site1 and Site2, respectively. The precession signals
for these frequencies should be readily observable with the present time
resolution of J-PARC limited by the Nyquist frequency $1/2\tau\simeq 6.25$ MHz
which is determined by the FWHM ($\tau\approx 80$ ns) of the pulsed $\mu^{+}$
beam. For comparison, an expected $\mu$SR spectrum at Site1 (0.84 MHz) is
shown in the inset of Fig. 1(a), where no such precession signal has been
observed.

Now, let us examine a more extreme possibility that Mu is delocalized along
the $c$-channel center due to the zero-point energy. For the estimation of
$E_{0}$ and $V_{\mathrm{b}}$, we investigated the total energy profile for $z$
position of H placed at the center of the $c$-channel in a $2\times 2\times 3$
superlattice, and evaluated its change $\Delta E_{\mathrm{tot}}$ with the
following equation,

$\displaystyle\Delta E_{\mathrm{tot}}(\bm{r})\equiv
E_{\mathrm{tot}}(\bm{r})-E_{\mathrm{min}},$ (3)

where $\bm{r}=(0,0.5,z)$ is the position vector in the unit cell and
$E_{\mathrm{min}}$ is the minimum value in $E_{\mathrm{tot}}(\bm{r})$. The $z$
dependence of $\Delta E_{\mathrm{tot}}(\bm{r})$ and the magnitude of $B_{\rm
loc}$ is shown in Fig. 4. $\Delta E_{\mathrm{tot}}(\bm{r})$ exhibits
oscillation with a $z/2$ period (shown on the left axis), with a maximum value
of $V_{\mathrm{b}}=0.353$ eV. In the harmonic approximation of the periodic
potential, the energy level splitting $\hbar\omega_{\mu}$ of Mu can be
determined from the relation
$\omega_{\mu}^{2}=2\pi^{2}V_{\mathrm{b}}/m_{\mu}d^{2}$ Yamaura  _et al._
(2019), where $m_{\mu}=105.658$ MeV/$c^{2}$ is the $\mu^{+}$ mass,
$d=c/2\simeq 0.155$ nm is the length corresponding to the potential period.
Then, $E_{0}=\tfrac{1}{2}\hbar\omega_{\mu}$ is estimated to be 0.163 eV, which
is about half of $V_{\mathrm{b}}$. Therefore, even when Mu is at the center of
the $c$-channel, it is expected to be localized at the minimum of $\Delta
E_{\mathrm{tot}}(\bm{r})$ corresponding to Site2 (the 4$c$ site). The $z$
dependence of the simulated $B_{\rm loc}$ for the reported AFM structure is
shown on the right axis of Fig. 4. We show the two cases where the Ru moment
$\lvert{\bm{m}}_{\mathrm{Ru}}\rvert=0.05\mu_{\mathrm{B}}$ with
${\bm{m}}_{\mathrm{Ru}}\parallel c$, and a 20∘ tilted to the (110) plane Zhu
_et al._ (2019). In the former case, $B_{\rm loc}$ cancels at $z=0.25$ (4$d$
site), but Mu is not likely to be stationary there because $\Delta
E_{\mathrm{tot}}(\bm{r})$ is at a maximum.

Figure 4: The total energy $\Delta E_{\mathrm{tot}}(\bm{r})$ and internal
magnetic field $B_{\rm loc}(\bm{r})$ along the $\bm{r}=(0,0.5,z)$ axis
(corresponding to the $c$-channel center) calculated in the $2\times 2\times
3$ superlattice of RuO2. The solid and dashed lines showing $B_{\rm
loc}(\bm{r})$ (mapped on the right axis) are for the cases where the Ru moment
${\bm{m}}_{\mathrm{Ru}}$ is either parallel to (001) (the $c$-axis) or it is
tilted 20∘ toward the (110) plane, respectively Zhu  _et al._ (2019). The
horizontal dashed line is the Mu zero-point energy.

Conversely, when the slow damping in the ZF spectrum at 5 K is interpreted as
part of a precession signal associated with extremely small $B_{\rm loc}$ that
appears with ordered Ru magnetic moments, we can estimate the maximum Ru
moment size from the precession frequency $f$ obtained by curve fit using the
following function,

$\displaystyle G_{z}(t)=A_{1}\left[\tfrac{1}{3}+\tfrac{2}{3}\cos(2\pi
ft+\varphi)\right],$ (4)

where $f=(\gamma_{\mu}/2\pi)B_{\mathrm{loc}}$, $\varphi$ (=0 in the present
experimental setup) is the initial phase of the precession,
$\gamma_{\mu}/2\pi=135.539$ MHz/T is the $\mu^{+}$ gyromagnetic ratio. Curve
fit of the ZF spectrum at 5 K shown in Fig. 1(a) yields $f=8.0(2)\times
10^{-3}$ MHz, corresponding to $B_{\mathrm{loc}}=5.9(2)\times 10^{-2}$ mT.
This value gives the upper boundary of $B_{\rm loc}$ for the field exerted
from Ru $4d$ electrons. By comparison with the lowest value of
$B_{\mathrm{loc}}$ among those obtained by the abovementioned simulation (5.96
mT at Site2), the upper limit for the $\lvert{\bm{m}}_{\mathrm{Ru}}\rvert$ can
be estimated to be $4.8(2)\times 10^{-4}\mu_{\mathrm{B}}$ for the reported
magnetic structures. Following the same procedure, the upper limit for
$\lvert{\bm{m}}_{\mathrm{Ru}}\rvert$ in the powder sample is also estimated to
be $6.1(3)\times 10^{-3}\mu_{\mathrm{B}}$ (see sm  for details). Value for
single crystal is the same order of magnitude as the nuclear magneton and is
nearly two orders of magnitude smaller than the value of
0.05$\mu_{\mathrm{B}}$ reported in the earlier study Berlijn  _et al._ (2017),
strongly disfavoring the presumed magnetic order in the bulk crystal. It
should be noted that these upper limits are even lower when the relaxation due
to the nuclear magnetic moments of Ru ($\approx 0.027$ $\mu$s-1 at the Site1)
is taken into account.

Finally, we briefly discuss a possibility of the slow relaxation due to the
fast fluctuating magnetic order. The fluctuation time $\tau_{c}$ can be
estimated using the Bloembergen-Purcell-Pound relation Bloembergen  _et al._
(1948), from the $\lambda_{0}\sim 0.007$ $\mu$s-1 in Eq. (2) and $B_{\rm
loc}\sim 6$ mT,

$\displaystyle\lambda_{0}(B_{\mathrm{LF}})\simeq
2(\gamma_{\mu}B_{\mathrm{loc}})^{2}\tau_{c}/\left(1+(\gamma_{\mu}B_{\mathrm{LF}}\tau_{c})^{2}\right),$

yielding $\tau_{c}\sim 0.13$ ns in ZF. While it appears static in neutron
scattering, it is detectable in $\mu$SR. However, $B_{\mathrm{LF}}$ required
to suppress the relaxation [$\lambda_{0}(B_{\mathrm{LF}})<0.001$ $\mu$s-1] is
estimated to be $>$20 T, which contradicts the experimental results, as shown
in Fig. 1(b). Thus, the possibility of the magnetic order with fast
fluctuation can be ruled out.

In summary, using $\mu$SR combined with DFT calculations, we provided the
experimental evidence that the electronic ground state of bulk crystal RuO2 is
a nonmagnetic metal. The $\mu$SR time spectrum at zero field shows only
monotonous relaxation without sinusoidal oscillation, indicating that no
homogeneous internal magnetic field is established at the Mu site(s). Our DFT
calculations ruled out the possibility for Mu to occupy sites where the
internal magnetic field accidentally cancels out for the reported AFM
structure. The upper limit of the Ru moment size is nearly two orders of
magnitude smaller than the reported value of 0.05$\mu_{\mathrm{B}}$. From
these results, we conclude that there is no bulk AFM order as reported in
RuO2.

###### Acknowledgements.

This work was supported by the MEXT Elements Strategy Initiative to Form Core
Research Center for Electron Materials (Grant No. JPMXP0112101001) and JSPS
KAKENHI (Grant No. 22K05275), and partially by the MEXT Program: Data Creation
and Utilization Type Material Research and Development Project under Grant No.
JPMXP1122683430. The authors would like to thank the help of the MUSE staff
during $\mu$SR experiment at MLF, J-PARC. The experiment was conducted under
the support of Inter-University Research Programs by Institute of Materials
Structure Science, KEK (Proposal No. 2019MS02).

## References

  * Over (2012) H. Over, Chem. Rev. 112, 3356 (2012), pMID: 22423981\.
  * Ryden and Lawson (2003) W. D. Ryden and A. W. Lawson, J. Chem. Phys. 52, 6058 (2003).
  * Sun  _et al._ (2017) Y. Sun, Y. Zhang, C.-X. Liu, C. Felser, and B. Yan, Phys. Rev. B 95, 235104 (2017).
  * Jovic  _et al._ (2018) V. Jovic, R. J. Koch, S. K. Panda, H. Berger, P. Bugnon, A. Magrez, K. E. Smith, S. Biermann, C. Jozwiak, A. Bostwick, E. Rotenberg, and S. Moser, Phys. Rev. B 98, 241101(R) (2018).
  * Berlijn  _et al._ (2017) T. Berlijn, P. C. Snijders, O. Delaire, H.-D. Zhou, T. A. Maier, H.-B. Cao, S.-X. Chi, M. Matsuda, Y. Wang, M. R. Koehler, P. R. C. Kent, and H. H. Weitering, Phys. Rev. Lett. 118, 077201 (2017).
  * Zhu  _et al._ (2019) Z. H. Zhu, J. Strempfer, R. R. Rao, C. A. Occhialini, J. Pelliciari, Y. Choi, T. Kawaguchi, H. You, J. F. Mitchell, Y. Shao-Horn, and R. Comin, Phys. Rev. Lett. 122, 017202 (2019).
  * Wadley  _et al._ (2016) P. Wadley, B. Howells, J. Železný, C. Andrews, V. Hills, R. P. Campion, V. Novák, K. Olejník, F. Maccherozzi, S. S. Dhesi, S. Y. Martin, T. Wagner, J. Wunderlich, F. Freimuth, Y. Mokrousov, J. Kuneš, J. S. Chauhan, M. J. Grzybowski, A. W. Rushforth, K. W. Edmonds, B. L. Gallagher, and T. Jungwirth, Science 351, 587 (2016).
  * Baltz  _et al._ (2018) V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Uchida  _et al._ (2020) M. Uchida, T. Nomoto, M. Musashi, R. Arita, and M. Kawasaki, Phys. Rev. Lett. 125, 147001 (2020).
  * Ruf  _et al._ (2021) J. P. Ruf, H. Paik, N. J. Schreiber, H. P. Nair, L. Miao, J. K. Kawasaki, J. N. Nelson, B. D. Faeth, Y. Lee, B. H. Goodge, B. Pamuk, C. J. Fennie, L. F. Kourkoutis, D. G. Schlom, and K. M. Shen, Nat. Commun. 12, 59 (2021).
  * $\mathrm{\check{S}}$mejkal  _et al._ (2020) L. $\mathrm{\check{S}}$mejkal, R. Gonz$\mathrm{\acute{a}}$lez-Hern$\mathrm{\acute{a}}$ndez, T. Jungwirth, and J. Sinova, Sci. Adv. 6, eaaz8809 (2020).
  * Feng  _et al._ (2022) Z. Feng, X. Zhou, L. Šmejkal, L. Wu, Z. Zhu, H. Guo, R. González-Hernández, X. Wang, H. Yan, P. Qin, X. Zhang, H. Wu, H. Chen, Z. Meng, L. Liu, Z. Xia, J. Sinova, T. Jungwirth, and Z. Liu, Nat. Electron. 5, 735 (2022).
  * González-Hernández  _et al._ (2021) R. González-Hernández, L. Šmejkal, K. Výborný, Y. Yahagi, J. Sinova, T. c. v. Jungwirth, and J. Železný, Phys. Rev. Lett. 126, 127701 (2021).
  * Bai  _et al._ (2022) H. Bai, L. Han, X. Y. Feng, Y. J. Zhou, R. X. Su, Q. Wang, L. Y. Liao, W. X. Zhu, X. Z. Chen, F. Pan, X. L. Fan, and C. Song, Phys. Rev. Lett. 128, 197202 (2022).
  * Bose  _et al._ (2022) A. Bose, N. J. Schreiber, R. Jain, D.-F. Shao, H. P. Nair, J. Sun, X. S. Zhang, D. A. Muller, E. Y. Tsymbal, D. G. Schlom, and D. C. Ralph, Nat. Electron. 5, 267 (2022).
  * Karube  _et al._ (2022) S. Karube, T. Tanaka, D. Sugawara, N. Kadoguchi, M. Kohda, and J. Nitta, Phys. Rev. Lett. 129, 137201 (2022).
  * Zhou  _et al._ (2021) X. Zhou, W. Feng, X. Yang, G.-Y. Guo, and Y. Yao, Phys. Rev. B 104, 024401 (2021).
  * Smolyanyuk  _et al._ (2023) A. Smolyanyuk, I. I. Mazin, L. Garcia-Gassull, and R. ValentÃ, arXiv (2023), arXiv:2310.06909 .
  * Nishida  _et al._ (1987) N. Nishida, H. Miyatake, D. Shimada, S. Okuma, M. Ishikawa, T. Takabatake, Y. Nakazawa, Y. Kuno, R. Keitel, J. H. Brewer, T. M. Riseman, D. L. Williams, Y. Watanabe, T. Yamazaki, K. Nishiyama, K. Nagamine, E. J. Ansaldo, and E. Torikai, Jpn. J. Appl. Phys. 26, L1856 (1987).
  * Uemura  _et al._ (1987) Y. J. Uemura, W. J. Kossler, X. H. Yu, J. R. Kempton, H. E. Schone, D. Opie, C. E. Stronach, D. C. Johnston, M. S. Alvarez, and D. P. Goshorn, Phys. Rev. Lett. 59, 1045 (1987).
  * Luke  _et al._ (1989) G. M. Luke, B. J. Sternlieb, Y. J. Uemura, J. H. Brewer, R. Kadono, R. F. Kiefl, S. R. Kreitzman, T. M. Riseman, J. Gopalakrishnan, A. W. Sleight, M. A. Subramanian, S. Uchida, H. Takagi, and Y. Tokura, Nature 338, 49 (1989).
  * Aeppli  _et al._ (1987) G. Aeppli, R. J. Cava, E. J. Ansaldo, J. H. Brewer, S. R. Kreitzman, G. M. Luke, D. R. Noakes, and R. F. Kiefl, Phys. Rev. B 35, 7129 (1987).
  * Uemura  _et al._ (1989) Y. J. Uemura, G. M. Luke, B. J. Sternlieb, J. H. Brewer, J. F. Carolan, W. N. Hardy, R. Kadono, J. R. Kempton, R. F. Kiefl, S. R. Kreitzman, P. Mulhern, T. M. Riseman, D. L. Williams, B. X. Yang, S. Uchida, H. Takagi, J. Gopalakrishnan, A. W. Sleight, M. A. Subramanian, C. L. Chien, M. Z. Cieplak, G. Xiao, V. Y. Lee, B. W. Statt, C. E. Stronach, W. J. Kossler, and X. H. Yu, Phys. Rev. Lett. 62, 2317 (1989).
  * Sonier  _et al._ (2000) J. E. Sonier, J. H. Brewer, and R. F. Kiefl, Rev. Mod. Phys. 72, 769 (2000).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Giannozzi  _et al._ (2020) P. Giannozzi, O. Baseggio, P. Bonf$\mathrm{\grave{a}}$, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, The Journal of Chemical Physics 152, 10.1063/5.0005082 (2020), 154105.
  * (28) See Supplemental Material at [URL will be inserted by publisher] for sample characterization, details of DFT calculations, simulations of internal magnetic field at Mu sites, and $\mu$SR results for a powder sample.
  * RodrÃguez-Carvajal (1993) J. RodrÃguez-Carvajal, Physica B: Condensed Matter 192, 55 (1993).
  * Dolomanov  _et al._ (2009) O. V. Dolomanov, L. J. Bourhis, R. J. Gildea, J. A. K. Howard, and H. Puschmann, J. Appl. Crystallogr. 42, 339 (2009).
  * Bourhis  _et al._ (2015) L. J. Bourhis, O. V. Dolomanov, R. J. Gildea, J. A. Howard, and H. Puschmann, Acta Crystallogr A Found Adv 71, 59 (2015).
  * Sheldrick (2015) G. M. Sheldrick, Acta Crystallographica Section C 71, 3 (2015).
  * (33) https://dalcorso.github.io/pslibrary.
  * Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Rogers  _et al._ (1969) D. B. Rogers, R. D. Shannon, A. W. Sleight, and J. L. Gillson, Inorg. Chem. 8, 841 (1969).
  * Kojima  _et al._ (2014) K. M. Kojima, T. Murakami, Y. Takahashi, H. Lee, S. Y. Suzuki, A. Koda, I. Yamauchi, M. Miyazaki, M. Hiraishi, H. Okabe, S. Takeshita, R. Kadono, T. Ito, W. Higemoto, S. Kanda, Y. Fukao, N. Saito, M. Saito, M. Ikeno, T. Uchida, and M. M. Tanaka, J. Phys.: Conf. Ser. 551, 012063 (2014).
  * Suter and Wojek (2012) A. Suter and B. Wojek, Phys. Proc. 30, 69 (2012).
  * (38) K. M. Kojima, private communication.
  * Uemura  _et al._ (1985) Y. J. Uemura, T. Yamazaki, D. R. Harshman, M. Senba, and E. J. Ansaldo, Phys. Rev. B 31, 546 (1985).
  * Hayano  _et al._ (1979) R. S. Hayano, Y. J. Uemura, J. Imazato, N. Nishida, T. Yamazaki, and R. Kubo, Phys. Rev. B 20, 850 (1979).
  * Cox  _et al._ (2006) S. F. J. Cox, J. L. Gavartin, J. S. Lord, S. P. Cottrell, J. M. Gil, H. V. Alberto, J. P. Duarte, R. C. Vil$\mathrm{\tilde{a}}$o, N. A. de Campos, D. J. Keeble, E. A. Davis, M. Charlton, and D. P. van der Werf, J. Phys.: Condens. Matter 18, 1079 (2006).
  * Ito  _et al._ (2023) T. U. Ito, W. Higemoto, and K. Shimomura, Phys. Rev. B 108, 224301 (2023).
  * Jadon  _et al._ (2021) A. Jadon, C. Rossi, M. Djafari-Rouhani, A. Est$\mathrm{\grave{e}}$ve, and D. Pech, Physics Open 7, 100059 (2021).
  * Marinopoulos  _et al._ (2018) A. G. Marinopoulos, R. C. Vil$\mathrm{\tilde{a}}$o, H. V. Alberto, and J. M. Gil, J. Phys.: Condens. Matter 30, 425503 (2018).
  * Hiraishi  _et al._ (2020) M. Hiraishi, K. M. Kojima, H. Okabe, S. Takeshita, A. Koda, R. Kadono, R. Khasanov, S. Iimura, S. Matsuishi, and H. Hosono, Phys. Rev. B 101, 174414 (2020).
  * Kim and Lai (2015) S. Kim and W. Lai, Comput. Mater. Sci. 102, 14 (2015).
  * Momma and Izumi (2011) K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * Yamaura  _et al._ (2019) J.-i. Yamaura, H. Hiraka, S. Iimura, Y. Muraba, J. Bang, K. Ikeuchi, M. Nakamura, Y. Inamura, T. Honda, M. Hiraishi, K. M. Kojima, R. Kadono, Y. Kuramoto, Y. Murakami, S. Matsuishi, and H. Hosono, Phys. Rev. B 99, 220505 (2019).
  * Bloembergen  _et al._ (1948) N. Bloembergen, E. M. Purcell, and R. V. Pound, Phys. Rev. 73, 679 (1948).


Si Lu and Ruisi Li 11institutetext: Portland State University

# DAC: Deep Autoencoder-based Clustering, a General Deep Learning Framework of
Representation Learning

Si Lu 11 Ruisi Li

###### Abstract

Clustering performs an essential role in many real world applications, such as
market research, pattern recognition, data analysis, and image processing.
However, due to the high dimensionality of the input feature values, the data
being fed to clustering algorithms usually contains noise and thus could lead
to in-accurate clustering results. While traditional dimension reduction and
feature selection algorithms could be used to address this problem, the simple
heuristic rules used in those algorithms are based on some particular
assumptions. When those assumptions does not hold, these algorithms then might
not work. In this paper, we propose DAC, Deep Autoencoder-based Clustering, a
generalized data-driven framework to learn clustering representations using
deep neuron networks. Experiment results show that our approach could
effectively boost performance of the K-Means clustering algorithm on a variety
types of datasets.

###### keywords:

clustering, K-Means, representation learning, deep neuron networks, deep
autoencoder

##  1 Introduction

Clustering is the task of grouping samples such that the ones in the same
group are more similar to each other than to the ones in other groups.
Nowadays, clustering performs as a basic and essential pre-processing step of
many real world applications. For example, it could be used to help with fake
news identification[6], document analysis[16], marketing and sales, etc.
Specifically, clustering algorithms can figure out useful information for the
applications via grouping according to a variety of data similarity metrics
and data grouping schemes. For example, similar patches could be used for
image denoising [3, 1, 2] or depth enhancement [9], and clustering could be
used to find good similar patches [8].

To let the samples be properly assigned to different groups(called clusters),
meaningful feature values of the samples need to be obtained first. However,
in real world applications, the data we get is often of high dimensions[5] and
usually contains noise, making the clustering difficult to succeed. For
example, in the MNIST dataset[7], each input hand-written digit image has 784
pixels. While we know some pixels (e.g. the ones at image corners) might not
be as useful as others(e.g. the ones around image centers), it is difficult to
manually distinguish them in clustering.

Traditional dimensionality reduction algorithms, namely Principle Component
Analysis(PCA)[10], Linear Discriminant Analysis(LDA)[4], and Canonical
Correlation Analysis(CCA)[13], could be used to reduce the number of features.
In addition, feature selection algorithms can be used to select from the
original feature values a set of useful and noiseless ones. These algorithms
aim to extract the core information given the redundant and correlated input
high-dimension data features. However, these algorithms often fail mainly due
to two reasons. Firstly, most of them require complex mathematical analysis,
which is difficult and time consuming as well. Secondly, their is no single
approach that could work for all types of datasets. Different datasets could
have different dimensions, data sizes and even might be used in totally
different applications. Some datasets are linear and some of them are non-
linear. As a result, it is difficult to find a way to generally work on all
types of datasets.

Recently, due to the emerging of the powerful deep neuron networks, deep
learning-based approaches have been introduced to learn better data
representations and achieve appealing performance improvements for clustering
algorithms. One simple approach is to learn representations using deep auto-
encoders. Specifically, the original input high dimensional features are fed
into a encoder that generates a low dimensional output. This output is further
fed to a decoder that tries to recover the raw input data as much as possible.
However, most of the existing approaches[11, 15] are using images as input and
thus using convolutional neuron networks in their work.

In this paper, we propose Deep Clustering Autoencoder, a simple but more
general framework for representation learning that takes feature vectors as
input. Thus, our approach could be applied to more generalized datasets. In
addition, according to the group labels, we propose a scheme to adaptively
weight all input features. We combine this estimated weight with the loss
function computation during training. Experiment results show that our
approach could effectively improve the performance of K-Means clustering
algorithm on different types of datasets, namely MNIST, Fashion-MNIST[14], as
well as Human Activities and Postural Transitions Data Set (HAPT)[12].

The rest of the paper is organized as follows: in section 2 we describe the
overview of our deep autoencoder-based clustering. We then describe the deep
autoencoder for representation learning in more details in section 3. We
finally show experimental results in section 4 and conclude in section 5.

Figure 1: Overview of our Deep Autoencoder-based Clustering on MNSIT dataset.
The autoencoder (consists of an encoder and a decoder) tries to encode and
decode the input features such that the decoded output is as close to the
input as possible. The input size is 28X28 = 784, the size of the learned low-
dimension representation is 10. In the testing stage, the learned encoder
output is then fed into the classic K-Means algorithm to do clustering.

##  2 Overview of Deep Autoencoder-based Clustering

Figure 1 shows an overview of our deep autoencoder-based clustering framework.
There are two main steps: training and clustering testing. In the training
step, a deep autoencoder with an encoder and a decoder is trained using the
training set. Here a flattened input vector is fed into the multilayer deep
encoder which has a low dimensional learned representation. This learned
representation is further fed into a decoder that tries to recover an output
of the same size as the input. The training process of this autoencoder tries
to reconstruct the input as much as possible. In the following clustering
step, we apply the autoencoder to the testing set. The output of the encoder
(learned representations) is then fed to a classic K-Means algorithm to do
clustering. The learned low dimensional representation vector contains key
information of the given input, and thus yield better clustering results.

##  3 Deep Autoencoder for Representation Learning

The architecture of our deep autoencoder for representation learning is shown
in Figure 1. As could be seen, the model is not as complex as some of the
advanced neuron networks. The reason is that we do not want our model to over-
fit in two-folds. First, we do not want our model to over-fit on the training
dataset over the testing dataset. Second, we do not want our model to over-fit
on the reconstruction problem it-self over the clustering problem. Thus, we
select a model of reasonable median complexity.

###  3.1 Encoder

The encoder aims to encode or compress the input data into a smaller size
representation, and at the same time preserve as much key information as
possible. As shown in Figure 1, the encoder consists of 8 layers, include the
input layer and the learned representation output layer. Here the input layer
is being normalized such that all its values is in the range of $(0,1)$.
Specifically, from the beginning, each larger layer is fully connected to the
next smaller layer followed by a couple of activation layers. There are mainly
two types of activation layers, Relu and Tanh, as shown in Equation 1 and 2.
Adding the Relu layers could introduce non-linearity to our model, making it
more robust against non-linear input data. The Tanh layer, on the other hand,
could transform the data into a normalized range of $(-1,1)$, to alleviate the
gradient vanishing/exploding problem.

$Relu(x)=max(0,x)$ (1)

$\begin{split}\cosh(x)&={\frac{{\rm e}^{x}+{\rm e}^{-x}}{2}}\\\
\sinh(x)&={\frac{{\rm e}^{x}-{\rm e}^{-x}}{2}}\\\
\tanh(x)&=\frac{\sinh(x)}{\cosh(x)}=\frac{{\rm e}^{x}-{\rm e}^{-x}}{{\rm
e}^{x}+{\rm e}^{-x}}\end{split}$ (2)

###  3.2 Decoder

The decoder aims to decode or decompress the encoded output to reconstruct the
original input data as much as possible. It contains nine layers, include the
input layer, which is the output of the encoder, and the final output layer.
Specifically, each smaller layer is fully connected to the next larger layer
followed by a Tanh activation layer. In addition, the decoder has a Sigmoid
activation layer (shown in Equation 3) at the final stage to enforce the
output values lie into the range of $(0,1)$.

$Sigmoid(x)=\frac{1}{1+e^{-x}}$ (3)

###  3.3 Objective Function

####  3.3.1 Clustering-weighted MSE Loss

While the goal of the classic autoencoder is to reconstruct the original input
as much as possible, it counts each input feature value equally. However, it
is possible that each individual input feature contributes differently to the
final clustering results. For example, in MNIST dataset, the pixels at the
four corners of almost all images are of the same color black (with zero
intensity input values), thus have no impact to the final clustering at all.
On the other hand, some pixels around the center of the images are likely to
perform more important roles. We thus propose a scheme to compute a
clustering-weighted MSE loss to let the autoencoder focus more on the
reconstruction of more important input feature values, as shown in Equation 4.

$L_{cmse}=\frac{\sum_{i=1}^{n}w_{i}(y_{i}-\hat{y_{i}})^{2}}{n}$ (4)

Here $w_{i}$ is the weight of each feature. It is computed using all ith
feature values sampled from a subset of the training dataset with $m$ samples.
Denote all ith feature values as $\\{x_{i}k|k=1,2,..,m\\}$ and the
corresponding ground truth group/cluster labels of the $m$ samples
$\\{l_{k}|k=1,2,..,m\\}$. The corresponding feature weight will be large if
both of the two following conditions are met. First, all sampled values in the
same groups/clusters have small differences. Second, all sampled values in
different groups/clusters have large differences. Thus, the weight is computed
as:

$w_{i}=\frac{\sum\limits_{l_{p}=l_{q}}e^{-(x_{ip}-x_{iq})^{2}}}{\sum\limits_{l_{p}=l_{q}}1}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\frac{\sum\limits_{l_{p}\neq
l_{q}}(1-e^{-(x_{ip}-x_{iq})^{2})}}{\sum\limits_{l_{p}\neq l_{q}}1}$ (5)

Figure 2: A map of the clustering weight computed for MNIST dataset using 1000
samples from the training set. It could be seen that pixels at boundaries and
corners are less important than the ones around image centers.

Figure 2 shows a map of the clustering weight computed for MNIST dataset using
1000 samples from the training set. Pixels at boundaries and corners are less
important than the ones around image centers, thus have smaller weights(white
means larger weights).

####  3.3.2 Final Objective Function

The final objective function then combines the Clustering-weighted MSE Loss
and a standard L2 norm regularization, as shown in Equation 6. Here the L2
norm regularization $L_{r}$ is computed using all parameters from the
autoencoder. $\beta$ is a balancing factor with a default value of $0.00001$.

$L=L_{cmse}+\beta\dot{L}_{r}$ (6)

##  4 Experimental Results

###  4.1 Dataset

We evaluate our approach on the classic MNIST hand-written digits dataset.
This dataset has $50,000$ images as the training set and $10,000$ images as
the testing set. There are 10 groups in total. We show some samples of MNIST
dataset in Figure 3.

Figure 3: Samples of the MNIST dataset.

###  4.2 Measurement Metrics

To evaluate our framework, we apply our trained encoder to the testing
dataset. We then compare the generated representations from our trained
encoder to the raw input features by applying them to the K-Means algorithm.
To measure the performance of clustering algorithms, we use the Adjusted Rand
Index (ARI). Specifically, this metrics computes a similarity between two
clustering results by considering all pairs of samples and counting pairs that
are assigned in the same or different clusters in the predicted and ground
truth clustering results. The proposed approach is denoted as DAC.

###  4.3 Experiment setup

We implement our framework in Python and PyTorch and test it on a desktop with
RTX 2080-Ti. We train the autoencoder for 200 epochs using Adam Optimization
Algorithm. The initial learning rate is set to 0.003 and will decrease with
the number of epochs during training.

###  4.4 Results on MNIST

Table 1 shows the quantitative performance of the proposed approach in terms
of $ARI$. Comparing to the raw K-Means algorithm, our approach ($DAC$) boosts
the K-Means algorithm’s performace from $0.3477$ to $0.6624$, which is a
$90.50\%$ boost. We also show some of the reconstructed results by our trained
autoencoder in Figure 4. It shows that our trained autoencoder can properly
reconstruct the raw input hand-written digits.

Table 1: Clustering results on MNIST testing dataset.  | K-Mesn | DAC
---|---|---
$ARI$ | 0.3477 | 0.6624
Figure 4: Sample results of our trained autoencoder on MNIST dataset. Top: raw
input images. Bottom: reconstructed images

###  4.5 Results on other datasets

To test the robustness of our approach against different data types, we apply
our method to two other datasets: Fashion-MNIST[14], and Human Activities and
Postural Transitions Data Set (HAPT)[12].

Fashion-MNIST is a similar dataset to MNIST, with the same image format and
image size. It has $60,000$ images as training set and $10,000$ images as
testing set. The only difference is the content: it contains images of 10
types of clothes. The ten categories are shown in Table 2. We show some
samples of this dataset in Figure 5.

Table 2: Fashion-MNIST category labels.  T-shirt/top |  Trouser |  Pullover |  Dress |  Coat
---|---|---|---|---
Sandal |  Shirt |  Sneaker |  Bag |  Ankle boot
Figure 5: Samples of the Fashion-MNIST dataset.

Human Activities and Postural Transitions Data Set is a dataset that has been
captured by smart phone’s sensors [12]. The authors captured 3-axial linear
acceleration and 3-axial angular velocity at a constant rate of 50Hz using the
embedded accelerometer and gyroscope of the device, which is a smartphone
(Samsung Galaxy S II). There are 30 volunteers whose ages are in the range of
19-48 years old. In their data capturing experiment, the volunteers was doing
one of twelve activities. There are six basic activities: three static
postures (standing, sitting, lying) and three dynamic activities (walking,
walking downstairs and walking upstairs). Another six postural transitions
that occurred between the static postures have also been added to the dataset.
These are: stand-to-sit, sit-to-stand, sit-to-lie, lie-to-sit, stand-to-lie,
and lie-to-stand. All twelve types of activities are shown in Table 3.

Table 3: HAPT category labels.  walking |  walking upstairs |  walking downstairs |  sitting
---|---|---|---
standing |  laying |  stand to sit |  sit to stand
standing |  laying |  stand to sit |  sit to stand
sit to lie |  lie to sit |  stand to lie |  lie to stand

The sensor signals (accelerometer and gyroscope) were then denoised by some
noise filters. The authors then sampled in fixed-width sliding windows of 2.56
sec and $50\%$ overlap (128 readings/window), leading to a sample size of 561
features. Each sample is captured when the volunteer is doing one type of
activities. During the capture process, $70\%$ of the volunteers were randomly
selected to generate the training set and $30\%$ were selected to generate the
testing set. In total, this dataset has 7767 samples for training and 3162
samples for testing.

Figure 6: Sample results of our trained on Fashion-MNIST dataset. Top: raw
input images. Bottom: reconstructed images

We apply our method to Fashion-MNIST dataset and report the results in Table
4. Here as the Fashion-MNIST is a more complex dataset, we modified the
autoencoder and show the modified autoencoder architecture in Figure 8. It can
be seen that comparing to using raw input features in K-Means clustering, our
method boosts ARI from 0.3039 to 0.4702, yields to a improvement of $54.7\%$.

We then apply our method to the HAPT dataset and report the results in Table
3. Here as this dataset’s inputs are of lower dimension than MNIST, we
modified the autoencoder accordingly and show the modified autoencoder
architecture in Figure 8. It can be seen that even with this temporal sequence
dataset, our method could effectively improve the K-Means algorithm’s
performance by $30\%$. These results also show that our method could be
generally applied to other data types. We also show some of the reconstructed
results by our trained autoencoder in Figure 6. It shows that our trained
autoencoder can properly reconstruct the raw input fashion images.

Table 4: Clustering results on Fashion-MNIST testing dataset.  | K-Mesn | DAC
---|---|---
$ARI$ | 0.3039 | 0.4702
Table 5: Clustering results on HAPT testing dataset.  | K-Mesn | DAC
---|---|---
$ARI$ | 0.4290 | 0.5594

##  5 Conclusion

In this paper, we propose DAC, Deep Autoencoder-based Clustering, a
generalized data-driven framework to learn low dimensional clustering
representations using trained deep neuron networks. Specifically, we train a
multi-layer deep autoencoder to encode and decode the raw input samples. The
encoded output of the encoder is then fed to a classic K-Means algorithm to do
clustering. We design a scheme to compute a clustering-based weight in the
training objective function to train the autoencoder and let it focus more on
the reconstruction of more important features. Experimental results show that
our approach could effectively boost the performance of a classic clustering
algorithm: K-Means by $30\%$ to $90\%$ on MNIST dataset. In addtion, our
method could be also applied to other types of clustering datasets, such as
Fashion-MNIST and Human Activities and Postural Transitions Data Set (HAPT).
Experimental results show that our framework could still be able to improve
K-Means algorithm’s performance by as much as $55\%$

Figure 7: Overview of our Deep Autoencoder-based Clustering on Fashion-MNSIT
dataset.

Figure 8: Overview of our Deep Autoencoder-based Clustering on HAPT dataset.

## References

  * [1] Buades, A., Coll, B., Morel, J.: A non-local algorithm for image denoising.  In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), vol. 2, pp. 60–65 (2005)
  * [2] Chen, F., Zhang, L., Yu, H.: External patch prior guided internal clustering for image denoising.  In: IEEE International Conference on Computer Vision (ICCV), pp. 603–611 (2015)
  * [3] Dabov, K., Foi, A., Katkovnik, V., Egiazarian, K.: Image denoising by sparse 3-d transform-domain collaborative filtering.  IEEE Transactions on Image Processing 16(8), 2080–2095 (2007)
  * [4] Duda, R.O., Hart, P.E., Stork, D.G.: Pattern classification. john wiley & sons.  Inc., New York 2 (2001)
  * [5] Han, J., Pei, J., Kamber, M.: Data mining: concepts and techniques.  Elsevier (2011)
  * [6] Hosseinimotlagh, S., Papalexakis, E.E.: Unsupervised content-based identification of fake news articles with tensor decomposition ensembles.  In: Proceedings of the Workshop on Misinformation and Misbehavior Mining on the Web (MIS2) (2018)
  * [7] LeCun, Y., Bottou, L., Bengio, Y., Haffner, P.: Gradient-based learning applied to document recognition.  Proceedings of the IEEE 86(11), 2278–2324 (1998)
  * [8] Lu, S.: Good similar patches for image denoising.  In: 2019 IEEE Winter Conference on Applications of Computer Vision (WACV), pp. 1886–1895. IEEE (2019)
  * [9] Lu, S., Ren, X., Liu, F.: Depth enhancement via low-rank matrix completion.  In: Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3390–3397 (2014)
  * [10] Pearson, K.: Liii. on lines and planes of closest fit to systems of points in space.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2(11), 559–572 (1901)
  * [11] Pu, Y., Gan, Z., Henao, R., Yuan, X., Li, C., Stevens, A., Carin, L.: Variational autoencoder for deep learning of images, labels and captions.  Advances in neural information processing systems 29, 2352–2360 (2016)
  * [12] Reyes-Ortiz, J.L., Oneto, L., Samà, A., Parra, X., Anguita, D.: Transition-aware human activity recognition using smartphones.  Neurocomputing 171, 754–767 (2016)
  * [13] Sun, Q.S., Zeng, S.G., Liu, Y., Heng, P.A., Xia, D.S.: A new method of feature fusion and its application in image recognition.  Pattern Recognition 38(12), 2437–2448 (2005)
  * [14] Xiao, H., Rasul, K., Vollgraf, R.: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  arXiv preprint arXiv:1708.07747 (2017)
  * [15] Yang, X., Deng, C., Zheng, F., Yan, J., Liu, W.: Deep spectral clustering using dual autoencoder network.  In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (2019)
  * [16] Zhao, Y., Karypis, G.: Evaluation of hierarchical clustering algorithms for document datasets.  In: Proceedings of the eleventh international conference on Information and knowledge management, pp. 515–524 (2002)


$\\{\epsilon\\}\cup\\{\langle{\varphi}\rangle+\bar{\varphi}\mid\varphi\in{L}|_{\not\Relbar\joinrel\mathrel{|}{s}},\bar{\varphi}\in\text{collect}_{\not\Relbar\joinrel\mathrel{|}}(M[\varphi],s)\\}$
collect⊑($(L,M),\bar{\psi}$) =
$\\{\epsilon\\}\cup\\{\langle{\varphi}\rangle+\bar{\varphi}\mid
i\in[{\bar{\psi}}],\varphi\in{L}|_{\sqsubseteq{\psi_{i}}},\bar{\varphi}\in\text{collect}_{\sqsubseteq}(M[\varphi],\langle{\psi_{j}\mid
j\neq i\wedge\psi_{j}\geq\varphi}\rangle_{j=1}^{{|{\bar{\psi}}|}})\\}$
${{\bm{\wedge}}_{\omega}[{\mathcal{L}}]}$ Fields:
$S:\text{Set}[{{{\bm{\wedge}}_{\omega}[{\mathcal{L}}]}}]$
$L:\text{LSet}[{\mathcal{L}}]$ Invariants: $\psi\in
L\Leftrightarrow\exists\bar{\varphi},i\in[{\bar{\varphi}}].\,(\bigwedge\bar{\varphi})\in
S\wedge\psi=\varphi_{i}$ Implementing
${}|_{\not\Relbar\joinrel\mathrel{|}{s}}$: $\\{\bigwedge\bar{\varphi}\in
S\mid\exists
i\in[{\bar{\varphi}}].\,\varphi_{i}\in{L}|_{\not\Relbar\joinrel\mathrel{|}{s}}\\}$
Implementing ${}|_{\sqsubseteq{\bigwedge\bar{\psi}}}$:
$\\{\bigwedge\bar{\varphi}\in S\mid\forall i\in[\bar{\psi}].\,\exists
j\in[{\bar{\varphi}}].\,\varphi_{j}\in{L}|_{\sqsubseteq{\psi_{i}}}\\}$
${{\bm{\exists}}_{X}[{\mathcal{L}}]}$ Fields:
$S:\text{Set}[{{{\bm{\exists}}_{X}[{\mathcal{L}}]}}]$
$L:\text{LSet}[{\mathcal{L}}]$ Invariants: $\varphi\in
L\Leftrightarrow(\exists X.\varphi)\in S$ Implementing
${}|_{\not\Relbar\joinrel\mathrel{|}{((\mathcal{U},\mathcal{I}),\mu)}}$:
$\\{\exists
X.\varphi\mid\varphi\in\bigcap_{\nu:X\to\mathcal{U}}{L}|_{\not\Relbar\joinrel\mathrel{|}{((\mathcal{U},\mathcal{I}),{\mu}\overleftarrow{\cup}{\nu})}}\\}$
Implementing ${}|_{\sqsubseteq{\exists X.\psi}}$:
$\bigcup_{\pi\in\mathfrak{S}_{X}}\\{\exists
X.\varphi\mid\varphi\in{L}|_{\sqsubseteq{\psi_{\pi}}}\\}$ where
$\psi_{\pi}=c_{\mathcal{L}}(\psi\pi)$ ${{\bm{\forall}}_{X}[{\mathcal{L}}]}$
Fields: $S:\text{Set}[{{{\bm{\forall}}_{X}[{\mathcal{L}}]}}]$
$L:\text{LSet}[{\mathcal{L}}]$ Invariants: $\varphi\in
L\Leftrightarrow(\forall X.\varphi)\in S$ Implementing
${}|_{\not\Relbar\joinrel\mathrel{|}{((\mathcal{U},\mathcal{I}),\mu)}}$:
$\\{\forall
X.\varphi\mid\varphi\in\bigcup_{\nu:X\to\mathcal{U}}{L}|_{\not\Relbar\joinrel\mathrel{|}{((\mathcal{U},\mathcal{I}),{\mu}\overleftarrow{\cup}{\nu}}}\\}$
Implementing ${}|_{\sqsubseteq{\forall X.\psi}}$:
$\bigcup_{\pi\in\mathfrak{S}_{X}}\\{\forall
X.\varphi\mid\varphi\in{L}|_{\sqsubseteq{\psi_{\pi}}}\\}$ where
$\psi_{\pi}=c_{\mathcal{L}}(\psi\pi)$ ${{\mathrlap{\bm{\exists}}\hskip
3.09998pt{\bm{\forall}}}_{X}[{\mathcal{L}}]}$ Fields:
$S:\text{Set}[{{{\mathrlap{\bm{\exists}}\hskip
3.09998pt{\bm{\forall}}}_{X}[{\mathcal{L}}]}}]$
$L_{\exists}:\text{LSet}[{\mathcal{L}}]$
$L_{\forall}:\text{LSet}[{\mathcal{L}}]$ Invariants: $\exists X.\varphi\in
L_{\exists}\Leftrightarrow(\exists X.\varphi)\in S$ $\forall X.\varphi\in
L_{\forall}\Leftrightarrow(\forall X.\varphi)\in S$ Implementing
${}|_{\not\Relbar\joinrel\mathrel{|}{s}}$:
${L_{\exists}}|_{\not\Relbar\joinrel\mathrel{|}{s}}\cup{L_{\forall}}|_{\not\Relbar\joinrel\mathrel{|}{s}}$
Implementing ${}|_{\sqsubseteq{QX.\psi}}$:
$\begin{cases}{L_{\forall}}|_{\sqsubseteq{\forall
X.\psi}}\cup{L_{\exists}}|_{\sqsubseteq{\exists X.\psi}}&\text{if
}Q=\exists\\\ {L_{\forall}}|_{\sqsubseteq{\forall X.\psi}}&\text{if
}Q=\forall\end{cases}$

We explain the implementation of $\text{LSet}[{\mathcal{L}}]$ assuming
standard data structures for sets and maps. In each case we construct the data
structure $\text{LSet}[{\mathcal{L}}]$ by augmenting a standard set data
structure, e.g., a hash-set, denoted $S:\text{Set}[{\mathcal{L}}]$, with
additional fields intended to facilitate the desired operations. Table 2 lists
the fields used in the $\text{LSet}[{\mathcal{L}}]$ data structure for each
language constructor, the invariants maintained, and the implementation of
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and
${\cdot}|_{\sqsubseteq{\psi}}$.

For the base case, $\text{LSet}[{\mathcal{L}_{A}}]$ is implemented without any
auxiliary fields, using straightforward iteration over the elements of $S$ for
computing ${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and
${\cdot}|_{\sqsubseteq{\psi}}$.

For ${{\bm{\vee}}[{\mathcal{L}_{1}},{\mathcal{L}_{2}}]}$, the data structure
has—in addition to $S$—two auxiliary data fields:
$L:\text{LSet}[{\mathcal{L}_{1}}]$, which holds the set of first-position
disjuncts present in $S$, and a map
$M:\text{Map}[{\mathcal{L}_{1}},{\text{LSet}[{\mathcal{L}_{2}}]}]$, which maps
each first-position disjunct to the set of its corresponding second-position
disjuncts (in our implementation $\text{Map}[{\cdot},{\cdot}]$ is a hash-map).
That is, the data structure maintains the invariants that a formula
$\varphi_{1}\vee\varphi_{2}$ is in the set $S$ iff $\varphi_{2}\in
M[\varphi_{1}]$, and that $L$ contains the same formulas as the domain (the
set of keys) of $M$, which we denote $\mathbf{dom}({M})$. For
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$, both disjuncts of a formula
must not be satisfied by $s$; therefore, this filter can be computed by using
${L}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ to get all first-position disjuncts
that are not satisfied by $s$, and for each $\varphi_{1}$ of them computing
${M[\varphi_{1}]}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ to get the relevant
second-position disjuncts. Similarly, to get all formulas subsuming
$\psi_{1}\vee\psi_{2}$, both disjuncts need to be subsumed; therefore,
${L}|_{\sqsubseteq{\psi_{1}}}$ is used to get all first-position disjuncts
that subsume $\psi_{1}$, and for each $\varphi_{1}$ of them
${M[\varphi_{1}]}|_{\sqsubseteq{\psi_{2}}}$ is used to get the relevant
second-position disjuncts subsuming $\psi_{2}$. Note that $L$ may not be an
antichain, even if $S$ is. For example, if
$S=\\{\varphi_{1}\vee\varphi_{2},\psi_{1}\vee\psi_{2}\\}$ with
$\varphi_{1}\sqsubseteq_{\mathcal{L}_{1}}\psi_{1}$ but
$\varphi_{2}\not\sqsubseteq_{\mathcal{L}_{2}}\psi_{2}$, then $S$ is an
antichain but $L$ (which holds $\\{\varphi_{1},\psi_{1}\\}$) is not.

The implementation for ${{\bm{\wedge}}[{\mathcal{L}_{1}},{\mathcal{L}_{2}}]}$
is similar to the binary disjunction case for
${\cdot}|_{\sqsubseteq{\psi_{1}\wedge\psi_{2}}}$, since the subsumption
relation is defined identically. However,
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ requires only one of the
conjuncts to be unsatisfied by $s$. We therefore hold two set-map pairs as
above: $L_{1}:\text{LSet}[{\mathcal{L}_{1}}]$,
$M_{1}:\text{Map}[{\mathcal{L}_{1}},{\text{LSet}[{\mathcal{L}_{2}}]}]$ and
$L_{2}:\text{LSet}[{\mathcal{L}_{2}}]$,
$M_{2}:\text{Map}[{\mathcal{L}_{2}},{\text{LSet}[{\mathcal{L}_{1}}]}]$. All
first-position conjuncts that do not satisfy $s$ can be retrieved using
${L_{1}}|_{\not\Relbar\joinrel\mathrel{|}{s}}$, and collecting for each such
$\varphi_{1}$ all $\varphi_{1}\wedge\varphi_{2}$ where $\varphi_{2}\in
M_{1}[\varphi_{1}]$ results exactly in all formulas in the set where the first
conjunct is not satisfied by $s$. In a symmetrical manner we find all formulas
in the set where the second conjunct is not satisfied by $s$. As mentioned
above, the ${\cdot}|_{\sqsubseteq{\psi_{1}\wedge\psi_{2}}}$ filter is
implemented identically to the
${{\bm{\vee}}[{\mathcal{L}_{1}},{\mathcal{L}_{2}}]}$ case using $L_{1}$,
$M_{1}$.

The case of $\text{LSet}[{{{\bm{\vee}}_{k}[{\mathcal{L}}]}}]$ generalizes the
case of ${{\bm{\vee}}[{\mathcal{L}_{1}},{\mathcal{L}_{2}}]}$, in the sense
that all disjuncts need to be unsatisfied by $s$ in
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$; and all disjuncts need to
subsume a (distinct) disjunct in ${\cdot}|_{\sqsubseteq{\bigvee\bar{\psi}}}$.
Therefore, these operations are implemented by computing them for all
disjuncts iteratively. We use a tree structure whose edges are labeled with
$\mathcal{L}$-formulas, and whose paths are understood as the sequences of
disjuncts as they appear in formulas $\bigvee\bar{\varphi}\in S$. Since we
consider canonical formulas, we are ensured that disjuncts are ordered
according to $\leq_{\mathcal{L}}$ along paths in the tree. Each node in the
tree has an $\text{LSet}[{\mathcal{L}}]$ that stores the
$\mathcal{L}$-formulas labeling its outgoing edges, which represent the next
possible disjuncts in the sequence of $\mathcal{L}$-formulas from the root to
the node. This generalizes the $L$ field in the binary case, which maintains
the possible first-position disjuncts, i.e., the next disjuncts w.r.t the
empty sequence. Using this construction,
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and
${\cdot}|_{\sqsubseteq{\bigvee\bar{\psi}}}$ are implemented by tree
traversals, with recursive applications of
${L}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and ${L}|_{\sqsubseteq{\psi_{i}}}$
for the sets $L:\text{LSet}[{\mathcal{L}}]$ in the nodes, which are used to
retrieve the next possible disjuncts not satified by $s$, or the next possible
disjuncts subsuming some disjunct in $\bigvee\bar{\psi}$, respectively. In the
case of ${\cdot}|_{\sqsubseteq{\bigvee\bar{\psi}}}$, since each resulting
disjunct needs to subsume a distinct formula in $\bar{\psi}$, we remove
$\psi_{i}$ from $\bar{\psi}$ during the recursive traversal if it caused us to
go down that edge. As an optimization, when going down an edge labeled
$\varphi$, we also remove from $\bar{\psi}$ all formulas
$\psi_{i}<_{\mathcal{L}}\varphi$, since such formulas cannot be subsumed by
formulas greater or equal to $\varphi$.

The case of $\text{LSet}[{{{\bm{\wedge}}_{\omega}[{\mathcal{L}}]}}]$ is
simpler. In addition to the standard
$S:\text{Set}[{{{\bm{\wedge}}_{\omega}[{\mathcal{L}}]}}]$, we maintain a field
$L:\text{LSet}[{\mathcal{L}}]$ aggregating all conjuncts from all formulas in
$S$. The formulas $\bigwedge\bar{\varphi}$ which are not satisfied by $s$ are
those with some conjunct in ${L}|_{\not\Relbar\joinrel\mathrel{|}{s}}$; and
the formulas $\bigwedge\bar{\varphi}$ subsuming $\bigwedge\bar{\psi}$ are
those where for all $i\in[{\bar{\psi}}]$ there is a conjunct of
$\bigwedge\bar{\varphi}$ in ${L}|_{\sqsubseteq{\psi_{i}}}$. These checks are
implemented in a straightforward way. An alternative approach is to store
conjunctions in a tree structure similar to the LSet implementation for
${{\bm{\vee}}_{k}[{\mathcal{L}}]}$, and perform
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and
${\cdot}|_{\sqsubseteq{\bigwedge\bar{\psi}}}$ using appropriate traversals of
the tree (analoguously to the way UBTrees [11] are used for finding both _sub_
sets and _super_ sets). However, observe that the savings in the
${{\bm{\wedge}}_{\omega}[{\mathcal{L}}]}$ case would be much more modest. For
example, the implementation of ${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$
for ${{\bm{\vee}}_{k}[{\mathcal{L}}]}$ uses a recursive traversal of the tree
that prunes away subtrees if the formula leading to them is satisfied by $s$
(since any disjunction containing that formula is satisfied by $s$). For
${{\bm{\wedge}}_{k}[{\mathcal{L}}]}$, discarding these paths is not possible,
since it suffices for _some_ conjunct to be violated by $s$ for the entire
conjunction to be unsatisfied by it. Thus,
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ would require recursively
traversing the entire tree, and our experience is that the overhead of
maintaining the tree structure makes this approach inefficient overall. We
note, however, that the implementation of LSet for unbounded conjunction is
the least optimized component of our approach, and might be improved further
in future work.

The cases of ${Q_{X}[{\mathcal{L}}]}$ for $Q\in\\{\exists,\forall\\}$ are each
implemented using an LSet $L$ over $\mathcal{L}$ to store the bodies of
formulas, which is used to perform
${\cdot}|_{\not\Relbar\joinrel\mathrel{|}{s}}$ and
${\cdot}|_{\sqsubseteq{QX.\psi}}$ by iterating over all assignments to the
variables, or over all permutations of variables, respectively, and using the
corresponding operations for $\text{LSet}[{\mathcal{L}}]$. The case of
${{\mathrlap{\bm{\exists}}\hskip 3.09998pt{\bm{\forall}}}_{X}[{\mathcal{L}}]}$
then uses these two data structures in a straightforward way.

Table 2 explained above relies heavily on maps that use formulas as keys.
However, computing hashes of formulas can be expensive, and therefore an
efficient implementation of $\text{LSet}[{\mathcal{L}}]$ requires an
additional optimization. Rather than using formulas as keys, in our
implementation each $\text{LSet}[{\mathcal{L}}]$ maintains a bidirectional
mapping between formulas and integer IDs, and uses the IDs whenever possible,
converting back to formulas only when producing the externally-facing output.
For example, the filters applied to $L$ in the implementation of
$\text{LSet}[{{{\bm{\vee}}[{\mathcal{L}_{1}},{\mathcal{L}_{2}}]}}]$ return a
set of IDs that are then used as keys for $M$, avoiding unnecessary hash
computations. (Other approaches, such as memoization of hashes or hash
consing, can alternatively be used.)

##  Appendix 0.D Language Optimizations

[29] uses several techniques for reducing the size of first-order languages
without losing precision. [29, Lemma 6], for example, implies the equivalence

$\displaystyle\bar{Q}\bar{x}.\,\varphi(\bar{x}_{\forall},\bar{x}_{\exists})\vee(\psi_{1}(\bar{x}_{\forall},\bar{x}_{\exists})\wedge\psi_{2}(\bar{x}_{\forall}))$
$\displaystyle\equiv\Big{[}\bar{Q}\bar{x}.\,\varphi(\bar{x}_{\forall},\bar{x}_{\exists})\vee\psi_{1}(\bar{x}_{\forall},\bar{x}_{\exists})\Big{]}\wedge\Big{[}\bar{Q}\bar{x}.\,\varphi(\bar{x}_{\forall},\bar{x}_{\exists})\vee\psi_{2}(\bar{x}_{\forall})\Big{]}.$

where $\bar{Q}\bar{x}$ is some quantifier prefix, $\bar{x}_{\forall}$ are the
universally quantified variables in $\bar{x}$ and $\bar{x}_{\exists}$ are the
the existentially quantified variables in $\bar{x}$. This allows [29] to
significantly reduce the number of formulas considered, since any formula that
can be equivalently decomposed in this way is semantically redundant
throughout the search for an inductive invariant.

However, as we have confirmed with the authors of [29], the proof of [29,
Lemma 6] contains an over-generalization, and the lemma in its general form
does not hold. For example, consider the formula

$\exists x.\forall y.\ (x\neq y)\vee(p(y)\wedge q(x))$

and the structure $(\mathcal{U},\mathcal{I})$ where $\mathcal{U}=\\{a,b\\}$
and $\mathcal{I}(p)=\\{a\\}$, $\mathcal{I}(q)=\\{b\\}$. It is easy to verify
that $(\mathcal{U},\mathcal{I})\not\models\exists x.\forall y.\ (x\neq
y)\vee(p(y)\wedge q(x))$, but

$\displaystyle(\mathcal{U},\mathcal{I})\models\exists x.\forall y.\ (x\neq
y)\vee p(y)\quad\text{and}\quad(\mathcal{U},\mathcal{I})\models\exists
x.\forall y.\ (x\neq y)\vee q(x).$

The problem with [29, Lemma 6] is that it might not work when some of the
universally quantified variables are quantified after existentially quantified
ones. Nonetheless, the reduction afforded by this kind of observation is
significant, so it is desirable to amend it, which we do as follows.

###### Lemma 1

For any quantification sequence $\forall\bar{x}\,\bar{Q}\bar{y}$ and formulas
$\varphi(\bar{x},\bar{y})$, $\psi_{1}(\bar{x},\bar{y})$, and
$\psi_{2}(\bar{x})$, the following semantic equivalence holds

$\displaystyle\forall\bar{x}.\bar{Q}\bar{y}.\
\varphi(\bar{x},\bar{y})\vee(\psi_{1}(\bar{x},\bar{y})\wedge\psi_{2}(\bar{x}))$
$\displaystyle\equiv\Big{[}\forall\bar{x}.\bar{Q}\bar{y}.\
\varphi(\bar{x},\bar{y})\vee\psi_{1}(\bar{x},\bar{y})\Big{]}\wedge\Big{[}\forall\bar{x}.\bar{Q}\bar{y}.\
\varphi(\bar{x},\bar{y})\vee\psi_{2}(\bar{x})\Big{]}.$

Before providing a proof for this lemma, we demonstrate how it allows us to
reduce the languages we are dealing with without losing expressive power. In
Sec. 6.2 we choose $k$-pDNF as our quantifier-free language, specified as

${{\bm{\vee}}[{{{\bm{\vee}}_{n}[{\cdot}]}},{{{\bm{\vee}}_{k-1}[{{{\bm{\wedge}}_{\omega}[{\cdot}]}}]}}]}.$

This permits us to select different base languages of first-order literals for
the $k$-pDNF clause (${{\bm{\vee}}_{n}[{\cdot}]}$) and for the $k$-pDNF cubes
(${{\bm{\vee}}_{k-1}[{{{\bm{\wedge}}_{\omega}[{\cdot}]}}]}$). Lemma 1 shows
that for a bounded language with a quantifier structure starting with
$\forall\bar{x}$, including literals containing only variables from $\bar{x}$
in the cubes does not add to the expressive power of $k$-pDNF—since the
decomposition described in Lemma 1 is able to substitute them with a
conjunction of $k$-pDNF formulas where cubes contain no such literals. Note,
however, that this decomposition may replace one $k$-pDNF formula with
$k$-pDNF formulas that have extra literals in their $k$-pDNF clause, so this
restriction may be a tradeoff between the parameter $n$ and the size of the
base language for the cubes. However, we find that in practice, using this
restriction captures the languages necessary to prove safety much more
precisely.

The proof for Lemma 1 is presented below.

######  Proof (Lemma 1)

We prove an equivalent second-order formulation of the claim, where
existentially quantified variables in $\bar{y}$ are replaced by existential
quantification over fresh function variables $\bar{f}$, similarly to the
process of Skolemization:

$\displaystyle\exists\bar{f}.\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\big{(}\psi_{1}(\bar{x},\bar{f}(\bar{x},\bar{z}))\wedge\psi_{2}(\bar{x})\big{)}$
$\displaystyle\equiv\Big{[}\exists\bar{f}.\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{1}(\bar{x},\bar{f}(\bar{x},\bar{z}))\Big{]}\wedge\Big{[}\exists\bar{f}.\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{2}(\bar{x})\Big{]}.$

Specifically, $\bar{z}$ consists of the universally quantified variables in
$\bar{y}$, and for all $i\in[{\bar{f}}]$, $f_{i}(\bar{v})$ replaces the $i$-th
existentially quantified variable in $\bar{y}$, where $\bar{v}$ are all the
universally quantified variables in $\bar{x},\bar{z}$ preceding that
existentially quantified variable in the quantifier prefix. To simplify the
notation we denote these function applications by $\bar{f}(\bar{x},\bar{z})$,
even though each of them may only depend on _some_ of the variables in
$\bar{z}$ (but note that all of them are dependent on $\bar{x}$).

Using distributivity, the left-hand side of the above is equivalent to

$\exists\bar{f}.\ \Big{[}\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{1}(\bar{x},\bar{f}(\bar{x},\bar{z}))\Big{]}\wedge\Big{[}\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{2}(\bar{x})\Big{]},$

which entails the right-hand side of the claim, since if there exists an
assignment to $\bar{f}$ satisfying both conjuncts, for each of the conjuncts
there exists an assignment satisfying it.

It is left to show entailment in the other direction. Let
$\sigma=(\mathcal{U},\mathcal{I})$ be a structure and $\mu_{1},\mu_{2}$ two
assignments to $\bar{f}$ such that

$\displaystyle(\sigma,\mu_{1})$
$\displaystyle\models\Big{[}\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{1}(\bar{x},\bar{f}(\bar{x},\bar{z}))\Big{]},\text{
and}$ $\displaystyle(\sigma,\mu_{2})$
$\displaystyle\models\Big{[}\forall\bar{x},\bar{z}.\
\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{2}(\bar{x})\Big{]}.$

Define the assignment $\mu$ to $\bar{f}$ which for all
$\bar{a}\in\mathcal{U}^{{|{\bar{x}}|}}$ and
$\bar{b}\in\mathcal{U}^{{|{\bar{z}}|}}$ satisfies

$\mu(\bar{f})(\bar{a},\bar{b})=\begin{cases}\mu_{1}(\bar{f})(\bar{a},\bar{b})&\text{if
$(\sigma,[\bar{x}\mapsto\bar{a}])\models\psi_{2}(\bar{x})$},\\\
\mu_{2}(\bar{f})(\bar{a},\bar{b})&\text{otherwise}.\end{cases}$

Note that this is only possible because, as we have explained, all function
applications in $\bar{f}(\bar{x},\bar{z})$ are dependent on $\bar{x}$. If that
were not the case, whether
$(\sigma,[\bar{x}\mapsto\bar{a}])\models\psi_{2}(\bar{x})$ could not be used
to determine $\mu(\bar{f})(\bar{a},\bar{b})$ in this way. Similarly, we cannot
condition the value of $\mu(\bar{f})(\bar{a},\bar{b})$ on $\bar{b}$ (other
than in the way $\bar{f}(\bar{x},\bar{z})$ depends on $\bar{z}$), since
functions may partially depend on it. In part, this is why this lemma cannot
be generalized to [29, Lemma 6].

The assignment $\mu$ shows that $\sigma$ satisfies the left-hand of the
equivalence, i.e.,

$(\sigma,\mu)\models\forall\bar{x},\bar{z}.\
\Big{[}\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{1}(\bar{x},\bar{f}(\bar{x},\bar{z}))\Big{]}\wedge\Big{[}\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))\vee\psi_{2}(\bar{x})\Big{]},$

because for all $\bar{a}\in\mathcal{U}^{{|{\bar{x}}|}}$ and
$\bar{b}\in\mathcal{U}^{{|{\bar{z}}|}}$, if
$(\sigma,\mu[\bar{x}\mapsto\bar{a},\bar{z}\mapsto\bar{b}])\models\psi_{2}(\bar{x})$
then the right-hand conjunct is satisfied, and the left-hand conjunct is
satisfied because

$\mu(\bar{f})(\bar{a},\bar{b})=\mu_{1}(\bar{f})(\bar{a},\bar{b}).$

Otherwise,
$(\sigma,\mu[\bar{x}\mapsto\bar{a},\bar{z}\mapsto\bar{b}])\not\models\psi_{2}(\bar{x})$,
which together with

$\mu(\bar{f})(\bar{a},\bar{b})=\mu_{2}(\bar{f})(\bar{a},\bar{b})$

implies
$(\sigma,\mu[\bar{x},\mapsto\bar{a},\bar{z}\mapsto\bar{b}])\models\varphi(\bar{x},\bar{f}(\bar{x},\bar{z}))$,
so both conjuncts are satisfied.

Thus, entailment holds in both directions, and the equivalence is established.

##  Appendix 0.E Least Fixpoint Computation via Symbolic Abstraction

Alg. 2 presents our implementation of the symbolic abstraction algorithm [28,
26] for computing the least fixpoint of the best abstract transformer of a
transition system $(\iota,\tau)$ in the abstract domain
${\mathcal{P}(\mathcal{L})}$ ordered by $\supseteq$.

##### Background and Notation.

The transition relation $\tau$ of a transition system over signature $\Sigma$
is a formula over signature $\Sigma\cup\Sigma^{\prime}$, where
$\Sigma^{\prime}$ is the primed copy of $\Sigma$ (in which each symbol is
replaced with a fresh primed duplicate). The set of transitions consists of
the pairs of states (structures over $\Sigma$) that jointly satisfy $\tau$.
More precisely, given
$(\mathcal{U},\mathcal{I}),(\mathcal{U}^{\prime},\mathcal{I}^{\prime})\in\mathbf{structs}[{\Sigma}]$
with $\mathcal{U}=\mathcal{U}^{\prime}$, we denote by
$\langle{(\mathcal{U},\mathcal{I}),(\mathcal{U}^{\prime},\mathcal{I}^{\prime})}\rangle\in\mathbf{structs}[{\Sigma\cup\Sigma^{\prime}}]$
the structure with the same universe $\mathcal{U}=\mathcal{U}^{\prime}$, where
$\mathcal{I}$ determines the interpretation of symbols from $\Sigma$ and
$\mathcal{I}^{\prime}$ determines the interpretation of symbols from
$\Sigma^{\prime}$. The set of transitions consists of
$\langle{s,s^{\prime}}\rangle$ such that
$\langle{s,s^{\prime}}\rangle\models\tau$. The semantics of the transition
system $(\iota,\tau)$ is described by a (concrete) transformer
$\mathcal{T}:{\mathcal{P}(\mathbb{S})}\to{\mathcal{P}(\mathbb{S})}$ defined by

$\mathcal{T}(S)=\\{s^{\prime}\in\mathbb{S}\mid
s^{\prime}\models\iota\lor\exists s\in S.\
\langle{s,s^{\prime}}\rangle\models\tau\\}.$

The least fixpoint of $\mathcal{T}$, denoted $\text{lfp}\mathcal{T}$, is the
set of reachable states of $(\iota,\tau)$. The abstract domain
${\mathcal{P}(\mathcal{L})}$ induces an abstract transformer
${\mathcal{T}}^{\sharp}:{\mathcal{P}(\mathcal{L})}\to{\mathcal{P}(\mathcal{L})}$
defined by ${\mathcal{T}}^{\sharp}=\alpha\circ\mathcal{T}\circ\gamma$, which
over-approximates the semantics of $(\iota,\tau)$ in the sense that
$\text{lfp}\mathcal{T}\subseteq\gamma(\text{lfp}{\mathcal{T}}^{\sharp})$.

Input: Transition system $(\iota,\tau)$, bounded first-order language
$\mathcal{L}$

Output: $R$ such that ${\uparrow}{R}=\text{lfp}{\mathcal{T}}^{\sharp}$ where
${\mathcal{T}}^{\sharp}$ is the best abstract transformer for $(\iota,\tau)$
in the abstract domain ${\mathcal{P}(\mathcal{L})}$

1

2$R$ : LSet[$\mathcal{L}$] := $\\{\bot_{\mathcal{L}}\\}$;

3

4while  _ $\exists s^{\prime}.\ s^{\prime}\models\iota\wedge\neg(\bigwedge R)$
or $\exists s,s^{\prime}.\ \langle{s,s^{\prime}}\rangle\models\bigwedge
R\wedge\tau\wedge\neg(\bigwedge R)^{\prime}$ _ do

/* loop invariant: ${\uparrow}{R}\supseteq\text{lfp}{\mathcal{T}}^{\sharp}$ */

weaken($R$, $s^{\prime}$) ;

// according to Alg. 1

5

6return _$R$_ ;

Algorithm 2 Least Fixpoint Computation of ${\mathcal{T}}^{\sharp}$

#####

Given $(\iota,\tau)$ and $\mathcal{L}$, Alg. 2 computes the representation of
$\text{lfp}{\mathcal{T}}^{\sharp}$. To do so, it stores and updates a set $R$,
which is the representation of the corresponding abstract element computed by
symbolic abstraction. Alg. 2 initializes the set $R$ to the representation of
$\mathcal{L}$ (the least element in the abstract domain
${\mathcal{P}(\mathcal{L})}$), which is $\\{\bot_{\mathcal{L}}\\}$. Note that
${\uparrow}{R}=\mathcal{L}$. In each iteration, a CTI $s^{\prime}$ is found
and $R$ is updated to a representation of
${\uparrow}{R}\sqcup\alpha(\\{s^{\prime}\\})$ by applying _weaken_ (Alg. 1).
The algorithm therefore maintains the invariant that in each iteration, $R$ is
the representation of $\alpha(S)$, where $S$ is the set of CTIs obtained so
far, i.e., ${\uparrow}{R}=\alpha(S)$.

A CTI is computed by finding a state $s^{\prime}$ such that either
$s^{\prime}\models\iota\wedge\neg(\bigwedge R)$, or there exists $s$ such that
$\langle{s,s^{\prime}}\rangle\models\bigwedge R\wedge\tau\wedge\neg(\bigwedge
R)^{\prime}$, where $(\bigwedge R)^{\prime}$ denotes the formula over
signature $\Sigma^{\prime}$ obtained from $(\bigwedge R)$ by replacing every
symbol from $\Sigma$ with its primed copy. Recall that a CTI should in fact be
found for the abstract element represented by $R$, but since $\bigwedge
R\equiv\bigwedge{\uparrow}{R}$, it suffices to compute a CTI based on the
representation $R$.

##  Appendix 0.F Optimizing the Search for CTIs

The symbolic abstraction algorithm described in App. 0.E and used in Sec. 6
requires performing SMT queries in order to check the inductiveness of a set
of formulas $R$. Since solvers often struggle when given even tens of formulas
(with quantifier alternations, but still in EPR), we employ several
optimizations. First, instead of directly checking the inductiveness of
$\bigwedge R$, we can separately check that each formula $\varphi\in R$ is
relatively inductive w.r.t $\bigwedge R$. Second, before checking the relative
inductiveness of a formula, we check if it is implied by formulas that were
already proven relatively inductive (implication queries are in general
cheaper than inductiveness queries since they involve only a single copy of
the signature). The above can also be parallelized across formulas. Third, we
use incremental SMT queries, as described in [14]. Lastly, in parallel to the
SMT queries, we attempt to find CTIs via simulations from previous CTIs.