# Correct and accurate polymorphic energy ordering of transition-metal
monoxides obtained from semilocal and onsite-hybrid exchange-correlation
approximations

Arghya Ghosh<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology, Hyderabad, India  Subrata Jana  Corresponding author:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Chemistry &
Biochemistry, The Ohio State University, Columbus, OH 43210, USA  Manish K
Niranjan<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Hyderabad, India  Fabien Tran  Institute of Materials Chemistry,
Vienna University of Technology, Getreidemarkt 9/165-TC, A-1060 Vienna,
Austria  David Wimberger  Institute of Materials Chemistry, Vienna University
of Technology, Getreidemarkt 9/165-TC, A-1060 Vienna, Austria  Peter Blaha
Institute of Materials Chemistry, Vienna University of Technology,
Getreidemarkt 9/165-TC, A-1060 Vienna, Austria  Lucian A. Constantin  Istituto
di Nanoscienze, Consiglio Nazionale delle Ricerche CNR-NANO, 41125 Modena,
Italy  Prasanjit Samal  School of Physical Sciences, National Institute of
Science Education and Research, HBNI, Bhubaneswar 752050, India

###### Abstract

The relative energetic stability of the structural phases of common
antiferromagnetic transition-metal oxides (MnO, FeO, CoO, and NiO) within the
semilocal and hybrid density functionals are fraught with difficulties. In
particular, MnO is known to be the most difficult case for almost all common
semilocal and hybrid density approximations. Here, we show that the meta-
generalized gradient approximation (meta-GGA) constructed from the cuspless
hydrogen model and Pauli kinetic energy density (MGGAC) can lead to the
correct ground state of MnO. The relative energy differences of zinc-blende
($zb$) and rock-salt ($rs$) structures as computed using MGGAC are found to be
in nice agreement with those obtained from high-level correlation methods like
the random phase approximation or quantum Monte Carlo techniques. Besides, we
have also applied the onsite hybrid functionals (closely related to DFT+$U$)
based on GGA and meta-GGA functionals, and it is shown that a relatively high
amount of Hartree-Fock exchange is necessary to obtain the correct ground-
state structure. Our present investigation suggests that the semilocal MGGAC
and onsite hybrids, both being computationally cheap, as methods of choice for
the calculation of the relative stability of antiferromagnetic transition-
metal oxides having potential applications in solid-state physics and
structural chemistry.

##  I Introduction

Following the pioneering works of Hohenberg, Kohn, and Sham Hohenberg and Kohn
(1964); Kohn and Sham (1965), density functional theory (DFT) has become a
highly successful and indispensable tool for studying the geometry and
electronic structure of solid-state and condensed matter systems Burke (2012);
Jones (2015); Becke (2014). Over the years, DFT-based first-principles studies
have greatly contributed to the progress in various fields relevant to
technology like nanoscience and molecular electronics. However, the accuracy
and reliability of DFT depend crucially on the various approximations for the
exchange-correlation (xc) energy functional Perdew and Schmidt (2001), which
includes all the many-body effects beyond the Hartree approximation. During
the last few decades, several accurate approximations for the semilocal xc
functionals, which are the computationally cheapest methods in DFT, have been
proposed. However, their application to transition-metal oxides (TMOs)
compounds having open $d$-shell still remains challenging Svane and Gunnarsson
(1990); Peng and Lany (2013); Schrön  _et al._ (2010); Schiller  _et al._
(2015); Peng and Perdew (2017)

Here in this paper, we revisit the relative accuracy of different levels of xc
methods to predict the ground-state properties of a prototypical open
$d$-shell TMO, namely MnO, which is potentially very interesting in industrial
applications, e.g., photoelectrochemical water splitting Kanan and Carter
(2012); Toroker and Carter (2013), solar energy conversion Peng and Lany
(2012), or magneto-piezoelectric effect Gopal  _et al._ (2004). A large number
of theoretical studies have been carried out for MnO phases Peng and Lany
(2013); Schrön  _et al._ (2010); Schiller  _et al._ (2015); Archer  _et al._
(2011); Franchini  _et al._ (2005); Peng and Perdew (2017) that have led to
the identification of different polymorphic phases Schrön  _et al._ (2010):
rock-salt ($rs$), zinc-blende ($zb$), and wurtzite ($wz$). Concerning
magnetism, MnO is antiferromagnetic (AF) with ferromagnetic planes stacked
along the [111] (AF2) and [001] (AF1) directions for the $rs$ and $zb$
structures, respectively. Among the phases $rs$-AF2 and $zb$-AF1, the first
one is found as the most stable one according to experiment Roth (1958) as
well as the random phase approximation (RPA) Peng and Lany (2013) and diffuse
Monte Carlo (DMC) methods Schiller  _et al._ (2015), which are high-level ab
initio methods. However, this is not the case with the common generalized
gradient approximations (GGA) and hybrid DFT xc methods that incorrectly
predict the energy of the $zb$-AF1 phase to be lower than that of $rs$-AF2
Peng and Lany (2013); Schrön  _et al._ (2010); Schiller  _et al._ (2015); Peng
and Perdew (2017). However, GGA+$U$ (GGA with a Hubbard $U$ correction)
predicts the correct ground state of MnO, although a large and unphysical $U$
value is required Kanan and Carter (2012).

One may note that in addition to the aforementioned works, a very large number
of other studies have investigated the ground-state or electronic properties
of the $rs$-AF2 phase of MnO. A certain number of DFT and beyond DFT methods
are used, and this includes many semilocal DFT methods (see Refs. Peng and
Perdew, 2017; Zhang  _et al._ , 2020; Sai Gautam and Carter, 2018; Tran  _et
al._ , 2018, 2020 for recent works), DFT+$U$ Anisimov  _et al._ (1991); Schrön
_et al._ (2012), various types of hybrids Tran  _et al._ (2006); Marsman  _et
al._ (2008); Liu  _et al._ (2019); Jana  _et al._ (2018a); Jana and Samal
(2019); Jana  _et al._ (2020a, b), the self-interaction corrected local
density approximation Svane and Gunnarsson (1990), the optimized effective
potential method Solovyev and Terakura (1998); Engel and Schmid (2009), model
Hamiltonian approach Archer  _et al._ (2011); Erten  _et al._ (2011); Hossain
_et al._ (2021); Hossain and Bhattacharjee (2020), the quasi-particle $GW$
method Massidda  _et al._ (1997); Kotani and van Schilfgaarde (2008); Rödl
_et al._ (2009), and dynamical mean field theory (DMFT) Kuneš  _et al._
(2008); Mandal  _et al._ (2019); Mukherjee  _et al._ (2014). Nevertheless, it
is important to mention that the most accurate of all these methods, namely,
RPA, DMC, $GW$, and DMFT are computationally much more expensive than DFT-
based methods.

Hence, from the point of view of efficiency, the preferred methods are the
semilocal xc functionals in DFT. In particular, we mention that
SCAN+rVV10+$U$, which consists of the strongly constrained and appropriately
normed (SCAN) meta-GGA Sun  _et al._ (2015) combined with the rVV10 van der
Waals (vdW) Sabatini  _et al._ (2013) functional and a Hubbard $U$ correction,
has been found to perform well in case of TMOs Peng and Perdew (2017). It is
clearly admitted that the meta-GGA functionals generally perform better than
the GGAs in describing solid-state properties Sun  _et al._ (2015); Jana  _et
al._ (2018b); Patra  _et al._ (2020a); Jana  _et al._ (2021a, 2020c, 2020d);
Ghosh  _et al._ (2021a). Concerning MnO, it has been shown that the
polymorphic structural energy difference computed using SCAN+rVV10+$U$, where
$U$ is determined from linear response theory, agrees very well with that
obtained from DMC values Peng and Perdew (2017). Besides the SCAN-based
methods, several other meta-GGA functionals are also proposed and tested for
solid-state properties with consistently improved accuracy Tran  _et al._
(2016); Sengupta  _et al._ (2018); Shahi  _et al._ (2018); Peng  _et al._
(2016); Mo  _et al._ (2017); Jana  _et al._ (2018b, c, 2019); Patra  _et al._
(2019a, b, 2021a); Jana  _et al._ (2021b); Ghosh  _et al._ (2021a). Within
those recent meta-GGAs, there are the functionals constructed from the
cuspless hydrogen model [(r)MGGAC] Patra  _et al._ (2019a); Jana  _et al._
(2021b) that show potential promising accuracy for different challenging
problems in solid-state physics Patra  _et al._ (2021a); Jana  _et al._
(2021b); Ghosh  _et al._ (2021a). It is also quite an efficient semilocal
functional that can predict band gaps of bulk and layered solids with
reasonable accuracy Patra  _et al._ (2020b, 2021b); Jana  _et al._ (2021b);
Tran  _et al._ (2021); Ghosh  _et al._ (2021a).

Inspired by the promising performance of the (r)MGGAC functionals for solids,
in the present paper we investigate the polymorphic energy ordering of MnO,
FeO, CoO, and NiO using these methods, with a particular focus on MnO. Our
results will also be compared with those obtained from other DFT methods and
higher-level quantum methods.

This paper is organized as follows. In Sec. II we first briefly describe the
methods that we used for the calculations. Then, in Sec. III we present and
discuss the results obtained for the polymorphic energy ordering of different
TMOs, and Sec. IV presents the conclusions.

##  II Methods of calculation

Belonging to the semilocal levels of approximation, the functionals PBE Perdew
_et al._ (1996), SCAN Sun  _et al._ (2015), r2SCAN Furness  _et al._ (2020),
MGGAC Patra  _et al._ (2019a), and rMGGAC Jana  _et al._ (2021b) are
considered for our calculations. The meta-GGA functionals (r2)SCAN and
(r)MGGAC, which depend on the density ($n=\sum_{i}|\psi_{i}|^{2}$), gradient
of the density ($\nabla n$), and KS kinetic-energy density (KED)
($\tau=\frac{1}{2}\sum_{i}|\nabla\psi_{i}|^{2}$), are implemented in the
generalized KS (gKS) Seidl  _et al._ (1996); Yang  _et al._ (2016); Perdew
_et al._ (2017) scheme. The semilocal functionals can be expressed as

$\displaystyle E_{\text{xc}}$ $\displaystyle=$
$\displaystyle\int\epsilon_{\text{xc}}(n,\nabla n,\tau)d^{3}r$ (1)
$\displaystyle=$
$\displaystyle\int\epsilon_{\text{x}}^{\text{LDA}}(n)F_{\text{xc}}(r_{s},s,\alpha^{\text{iso}})d^{3}r,$

where
$\epsilon_{\text{x}}^{\text{LDA}}=-\left(3/4\right)\left(3/\pi\right)^{1/3}n^{4/3}$
is the exchange energy density of the local density approximation (LDA),
$F_{\text{xc}}$ is the xc enhancement factor (with no dependency on $\tau$ for
GGAs), and $r_{s}=\left(3/\left(4\pi n\right)\right)^{1/3}$ is the Wigner-
Seitz radius. Usually, the $\nabla n$\- and $\tau$-dependencies of
$F_{\text{xc}}$ are expressed via dimensionless variables, like the reduced
density gradient $s=|\nabla
n|/\left(2\left(3\pi^{2}\right)^{1/3}n^{4/3}\right)$ and the iso-orbital
indicator $\alpha^{\text{iso}}=(\tau-\tau^{\text{W}})/\tau^{\text{UEG}}$,
where $\tau^{\text{W}}$ is the von Weizsäcker KED and $\tau^{\text{UEG}}$ is
the KED of the uniform electron gas (UEG). Actually, it is important to note
that $\alpha^{\text{iso}}$ recognizes regions with single bonds, overlapping
orbitals, and uniform density Sun  _et al._ (2013); Della Sala  _et al._
(2016).

Figure 1: Enhancement factors $F_{\text{xc}}$ plotted as a function of $s$
(upper panel) and $\alpha^{\text{iso}}$ (lower panel). The value of the two
other variables (that are kept fixed) are indicated in the respective panels.
Note the different scales on the vertical axis.

Note that the SCAN functional was constructed such that it satisfies 17 exact
mathematical constraints that can be satisfied by a meta-GGA. For instance,
the exchange component of the xc enhancement factor, $F_{x}^{\text{SCAN}}$,
recovers the exact fourth-order gradient approximation of exchange. r2SCAN is
a slightly modified version of SCAN that alleviates numerical problems
encountered with SCAN. The development of the MGGAC functional is rather
unusual and unconventional. The exchange part is based on the Becke-Roussel
approachBecke and Roussel (1989), and with a cuspless hydrogen exchange hole
density Patra  _et al._ (2019c). MGGAC differs from SCAN for the following two
reasons. First, only the exchange component of MGGAC is a meta-GGA, while the
correlation part is a GGA. Second, $F_{x}^{\text{MGGAC}}$ depends only on
$\alpha^{\text{iso}}$ (and not on $s$). However, the correlation part of
rMGGAC is of the meta-GGA type Jana  _et al._ (2021b). It is worth to mention
that (r2)SCAN and (r)MGGAC respect the strongly tightened bound exchange
($F_{x}\leq 1.174$ Perdew  _et al._ (2014); Sun  _et al._ (2015)) and possess
ultranonlocality effects, which is important for the band gap problem
Aschebrock and Kümmel (2019); Patra  _et al._ (2020b, 2019a). It may also be
noted that MGGAC correlation is not free from the one-electron self-
interaction error, whereas rMGGAC is.

In Fig. 1, the enhancement factor $F_{\text{xc}}$ is shown for the PBE,
(r2)SCAN, and (r)MGGAC functionals. They are plotted as functions of $s$ or
$\alpha^{\text{iso}}$, with fixed chosen values for the other parameters. An
obvious difference between the GGA PBE and the meta-GGAs concerns the sign of
$\partial F_{\text{xc}}/\partial s$. While $\partial F_{\text{xc}}/\partial s$
is positive for PBE, it is negative for all four meta-GGAs. Note that in
general the slope $\partial F_{\text{xc}}/\partial s$ of meta-GGAs can be
positive or negative depending on the particular meta-GGA, but also, to a
lesser extent, on the chosen fixed values of $r_{s}$ or $\alpha^{\text{iso}}$
(see Ref. Tran  _et al._ , 2020 for plots of $F_{\text{xc}}$ for other meta-
GGAs). The other difference between PBE and the meta-GGAs concerns of course
the variation with respect to $\alpha^{\text{iso}}$. Since PBE is a GGA,
$\partial F_{\text{xc}}/\partial\alpha^{\text{iso}}=0$. The meta-GGAs have a
negative value of $\partial F_{\text{xc}}/\partial\alpha^{\text{iso}}$ as for
most meta-GGAs Tran  _et al._ (2020). As discussed in Ref. Aschebrock and
Kümmel, 2019, a more negative slope $\partial
F_{\text{xc}}/\partial\alpha^{\text{iso}}$ leads to a larger derivative
discontinuity, and consequently also to a larger band gap.

In addition to the aforementioned functionals, the onsite hybrid functionals,
which are proposed specifically for systems with strongly correlated electrons
Tran  _et al._ (2006), are also considered in this work. They share close
similarities with the DFT+$U$ method, in particular since they are also
applied only to the strongly correlated electrons, while the rest of the
electrons are treated at the semilocal level. However, there are a few
fundamental and technical differences. For instance, whereas $U$, which
represents the screened Coulomb interaction, is the parameter on which the
DFT+$U$ results depend, in the case of the onsite hybrids it is the amount of
Hartree-Fock exchange that can be tuned to vary the results. Also, in the
onsite hybrids the double counting correction is calculated at the semilocal
DFT level, while a more ad hoc correction is used in DFT+$U$. The two onsite
hybrid functionals that are tested in the present work, PBE-$\alpha$ and
SCAN-$\alpha$, use PBE and SCAN as underlying semilocal functionals,
respectively. For both, three different values of the amount of Hartree-Fock
exchange $\alpha$ (0.25, 0.40, and 0.55) are used. Technical details about the
onsite hybrids and their implementations can be found in Refs. Novák  _et al._
, 2006; Tran  _et al._ , 2006. Note that in the SCAN-$\alpha$ calculations,
the SCAN functional was applied non-self-consistently and the PBE potential
was used instead.

##  III Results

Table 1: Total energy difference $\Delta
E=E_{zb}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ (in meV/f.u.) between the $zb$-AF1
and $rs$-AF2 phases of MnO. A positive energy difference indicates that
$rs$-AF2 is energetically more stable than $zb$-AF1 (as experimentally
determined). The MGGAC, rMGGAC, r2SCAN, and onsite hybrids results are
calculated for this work. The results obtained with the other methods are from
Refs. Schiller  _et al._ , 2015; Peng and Lany, 2013; Peng and Perdew, 2017.
$U$ and $V$ values are in eV.

Functional type | Functional | $\Delta E$
---|---|---
Semilocal | PBE | -244a
| SCAN | -79a
| r2SCAN | -33
| MGGAC | 118
| rMGGAC | 101
Semilocal+vdW | PBE+TS | -151a
Semilocal(+vdW)+$U$ | PBE+$U$ ($U=3$) | -19b
| PBE+TS+$U$ ($U=2$) | 88b
| SCAN+rVV10+$U$ ($U=2.8$) | 138a
Hybrid | HSE06 | -28a
Onsite hybrids | PBE-$\alpha$ ($\alpha=0.25$) | -84
| PBE-$\alpha$ ($\alpha=0.40$) | -22
| PBE-$\alpha$ ($\alpha=0.55$) | 23
| SCAN-$\alpha$ ($\alpha=0.25$) | 38
| SCAN-$\alpha$ ($\alpha=0.40$) | 98
| SCAN-$\alpha$ ($\alpha=0.55$) | 145
Higher level methods | PBE+$V{{}_{GW}}$ ($U=7$, $V=3$) | 64b
| RPA@PBE+$U$ ($U=3$) | 67b
| RPA@PBE+$V{{}_{GW}}$ ($U=7$, $V=3$) | 131b
| DMC | 132c

33footnotetext: Ref. Schiller  _et al._ , 2015.

Figure 2: Total energy difference $\Delta
E=E_{zb}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ (in meV/f.u.) between the $zb$-AF1
and $rs$-AF2 phases of MnO. A positive energy difference indicates that
$rs$-AF2 is energetically more stable than $zb$-AF1 (as experimentally
determined). The MGGAC, rMGGAC, r2SCAN, and onsite hybrids results are
calculated for this work. The results obtained with the other methods are from
Refs. Schiller  _et al._ , 2015; Peng and Lany, 2013; Peng and Perdew, 2017.
Figure 3: Total energy versus volume (in Å3 per atom) for the $rs$-AF2 and
$zb$-AF1 phases of MnO obtained from different methods. The total energy is
set to zero at the minimum of the curve of the $rs$-AF2 phase.

Figure 4: Total energy difference $\Delta
E=E_{zb}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ (in meV/f.u.) between the $zb$-AF1
and $rs$-AF2 phases of FeO (upper panel), CoO (middle panel), and NiO (lower
panel). A positive energy difference indicates that $rs$-AF2 is energetically
more stable than $zb$-AF1. The MGGAC, rMGGAC, r2SCAN, and onsite hybrids
results are calculated for this work. See Refs. Schiller  _et al._ , 2015;
Peng and Lany, 2013; Peng and Perdew, 2017 for the detailed results obtained
with the other methods.

The calculated total energy difference $\Delta
E=E_{zb}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ and stability of the $zb$-AF1 phase
relative to the experimentally determined ground-state $rs$-AF2 phase of MnO
are shown in Table 1 and Fig. 2. The results obtained with PBE, (r2)SCAN,
(r)MGGAC, and the onsite hybrids are obtained for the present work and are
compared to results from the literature obtained with other methods. As
mentioned in Sec. I, methods from various levels of theory, including
semilocal, vdW-corrected, and high-level correlation methods, have been used
and proposed in previous works to correctly describe the relative phase
stability of MnO and to provide estimates of $\Delta E$. In the present work,
we consider the DMC value $\Delta E^{\text{DMC}}=132$ meV/f.u. from Ref.
Schiller  _et al._ , 2015 as the reference benchmark.

We mention again that RPA correctly leads to $rs$-AF2 as the ground-state
phase Peng and Lany (2013). Compared to DMC, the value $\Delta E=67$ meV/f.u.
from RPA@PBE+$U$ ($U=3$ eV) (PBE+$U$ orbitals used as input to RPA) is too
small by a factor of 2, while $\Delta E=131$ meV/f.u from RPA@PBE+$V_{GW}$
($U=7$ eV and $V=3$ eV) matches perfectly the DMC value (in PBE+$V_{GW}$ a
nonlocal external potential $V_{GW}$ on Mn $d$ orbitals is used on top of the
onsite $U$). In Ref. Peng and Perdew, 2017 it is shown that the correct phase
ordering can also be obtained by adding vdW (Tkatchenko-Scheffler (TS)
Tkatchenko and Scheffler (2009) or rVV10 Vydrov and Van Voorhis (2010)) and
Hubbard $U$ (obtained from linear-response approach) corrections to a
semilocal functional. The proposed PBE+TS+$U$ ($U=3.2$ eV) Peng and Perdew
(2017) and SCAN+$r$VV10+$U$ ($U=2.8$ eV) Peng and Perdew (2017) lead to
$\Delta E=88$ and 138 meV/f.u., respectively, the latter value agreeing very
well with DMC.

However, as clearly visible in Fig. 2, there are also a certain number of
other popular methods that fail in predicting the correct ground-state of MnO.
Besides plain PBE, it is also the case with PBE+$U$ with a small value of $U$
(3 eV), the hybrid functional HSE06, and PBE+TS Peng and Lany (2013). The
meta-GGAs (r2)SCAN also lead to the wrong phase for the ground state, however
the two other meta-GGAs tested in this work, MGGAC and rMGGAC, lead to
positive energy differences, indicating that $rs$-AF2 is more stable than
$zb$-AF1. Thus, MGGAC and rMGGAC are successful in predicting correctly the
polymorphic energy ordering of the MnO phases. These are very interesting
results, in particular when considering that no Hubbard $U$ or vdW corrections
are added to (r)MGGAC. Actually, as evident from the results, MGGAC and rMGGAC
give the best performance among the semilocal methods. They lead to values for
$\Delta E$ of $118$ and $101$ meV/f.u., respectively, which agree quite well
with the reference value 132 meV/f.u. from DMC.

Regarding the performance of the onsite hybrids PBE-$\alpha$ and
SCAN-$\alpha$, also tested for the present work, we can see that the energetic
ordering of the $rs$-AF2 and $zb$-AF1 phases depends on the amount of the
Hartree-Fock exchange $\alpha$. For a correct energetic ordering, PBE-$\alpha$
requires a value of $\alpha$ larger than 0.40, which is quite high and clearly
larger than the standard value $\alpha=0.25$. For SCAN-$\alpha$, $\alpha=0.25$
is already sufficient. The observed trend with the onsite hybrids is that the
more $\alpha$ is large, the more the energetic ordering will go in the right
direction. This is in agreement with the results from Schiller et al. Schiller
_et al._ (2015) who considered hybrid functionals based on PBE (however, they
found that $rs$-AF2 becomes more stable already at $\alpha=0.10$). This
behavior is also similar to what is obtained with DFT+$U$, with $U$ playing
the same role as $\alpha$ Kanan and Carter (2012); Peng and Lany (2013).

In Fig. 3, we show the total energy versus volume for the two phases of MnO
obtained with various methods. As discussed above, unlike for the PBE and SCAN
methods the total-energy volume curve of the $rs$-AF2 phase is below that for
the $zb$-AF1 phase when MGGAC and rMGGAC are used.

Table 2: Total energy difference $\Delta E=E_{zb}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ (in meV/f.u.) between the $zb$-AF1 and $rs$-AF2 phases of MnO obtained non-self-consistently with different functionals using the PBE orbitals and electron density. The calculations are performed using the all-electron code WIEN2k. | PBE | SCAN@PBE | r2SCAN@PBE | MGGAC@PBE | rMGGAC@PBE
---|---|---|---|---|---
$\Delta E$ | -243 | -87 | -57 | 109 | 84
Table 3: Lattice constant $a_{0}$, bulk modulus $B_{0}$, band gap $E_{g}$, and
magnetic moment $\mu$ of the Mn atom of the $rs$-AF2 phase of MnO obtained
from different methods. The calculations with the semilocal functionals are
performed using the VASP code, while the onsite hybrid calculations are done
using the WIEN2k code.

Method |  $a_{0}$ (Å) |  $B_{0}$ (GPa) |  $E_{g}$ (eV) |  $\mu$ ($\mu_{\text{B}}$)
---|---|---|---|---
PBE | 4.438 | 149 | 0.72 | 4.38
SCAN | 4.411 | 163 | 1.47 | 4.49
r2SCAN | 4.418 | 166 | 1.52 | 4.50
MGGAC | 4.381 | 185 | 1.77 | 4.52
rMGGAC | 4.392 | 178 | 1.70 | 4.52
PBE-$\alpha$ ($\alpha=0.25$) | 4.488 | 144 | 1.21 | 4.65
SCAN-$\alpha$ ($\alpha=0.25$)111Ref. Peng and Perdew, 2017. | 4.447 | 161 | N/A | N/A
PBE-$\alpha$ ($\alpha=0.40$) | 4.512 | 141 | 1.34 | 4.74
SCAN-$\alpha$ ($\alpha=0.40$)111Since the SCAN functional in SCAN-$\alpha$ was applied non-self-consistently (the PBE potential was used) the values of $E_{g}$ and $\mu$ are omitted. | 4.465 | 158 | N/A | N/A
PBE-$\alpha$ ($\alpha=0.55$) | 4.531 | 139 | 1.45 | 4.80
SCAN-$\alpha$ ($\alpha=0.55$)111Since the SCAN functional in SCAN-$\alpha$ was applied non-self-consistently (the PBE potential was used) the values of $E_{g}$ and $\mu$ are omitted. | 4.482 | 156 | N/A | N/A
Expt. 222Ref. Peng and Lany, 2013. | 4.4365, 4.4315, 4.4302 | 149.6, 146.7, 148/144 | 3.6-3.9 | 4.58

22footnotetext: See Ref. Schrön  _et al._ , 2010 and references therein.

Table 4: Values of $E_{\text{xc}}$ in the $rs$-AF2 and $zb$-AF1 phases of MnO calculated with the SCAN and MGGAC functionals. The results are obtained non-self-consistently (with the PBE electron density) at the SCAN equilibrium geometry using the WIEN2k code. The total value in the cell is decomposed into the Mn and O atomic spheres and interstitial. The differences between the two phases are also shown, as well as their difference between the two functionals (last column). The units are in Ry/f.u. for columns 1, 2, 4, and 5 and in meV/f.u. for columns 3, 6, and 7. | SCAN | MGGAC | SCAN$-$MGGAC
---|---|---|---
| $rs$ | $zb$ | $rs-zb$ | $rs$ | $zb$ | $rs-zb$ | $rs-zb$
Cell | -122.693 | -122.595 | -1324 | -123.077 | -122.966 | -1515 | 191
Mn | -104.304 | -104.246 | -786 | -104.470 | -104.409 | -829 | 43
O | -16.887 | -16.934 | 646 | -17.105 | -17.152 | 649 | -3
Interstitial | -1.502 | -1.415 | -1183 | -1.502 | -1.404 | -1335 | 152
Figure 5: Two-dimensional plot in MnO of the difference
$\epsilon_{\text{xc}}^{\text{SCAN}}-\epsilon_{\text{xc}}^{\text{MGGAC}}$
between the SCAN and MGGAC xc energy density [see Eq. (1)] obtained with
WIEN2k from non-self-consistent calculations using the PBE electron density.
The left and right panels show the (110) plane of $rs$-AF2 and the (001) plane
of $zb$-AF1, respectively. The blue and red regions correspond to negative and
positive values, respectively. Figure 6: Two-dimensional plot in MnO of the
difference $F_{\text{xc}}^{\text{SCAN}}-F_{\text{xc}}^{\text{MGGAC}}$ between
the SCAN and MGGAC xc enhancement factors [see Eq. (1)] obtained with WIEN2k
from non-self-consistent calculations using the PBE electron density. The left
and right panels show the (110) plane of $rs$-AF2 and the (001) plane of
$zb$-AF1, respectively. The blue and red regions correspond to negative and
positive values, respectively. The atoms are the same as indicated in Fig. 5.

Figure 7: Two-dimensional plots in MnO of $r_{s}$ (upper panels), the
$\sigma$-spin reduced density gradient $s_{\sigma}$ (middle panels), and the
$\sigma$-spin iso-orbital indicator $\alpha_{\sigma}^{\text{iso}}$ (lower
panels). The left and right panels show the (110) plane of $rs$-AF2 and the
(001) plane of $zb$-AF1, respectively. The atoms are the same as indicated in
Fig. 5.

The (r2)SCAN and (r)MGGAC results discussed so far are obtained self-
consistently using the Vienna Ab initio Simulation Package (VASP) code Kresse
and Hafner (1993); Kresse and Furthmüller (1996, 1996); Kresse and Joubert
(1999). It may be of interest to determine how important is self-consistency
for the results, in particular in the case of (r)MGGAC, which lead to very
good results. It is generally believed that in most cases self-consistency
plays only a minor role and that using for instance the PBE electron density
and orbitals (instead of those from (r)MGGAC) would not affect much the total
energy computed with the (r)MGGAC energy functional. However, it may not be
always the case, as discussed in Ref. Kim  _et al._ , 2013 for the closely
related density-driven error. To make this point more clear for the results in
the present work, we also calculated $\Delta E$ for (r)MGGAC and (r2)SCAN
using the PBE orbitals and electron density. These non-self consistent
calculations are performed using the all-electron code WIEN2kBlaha  _et al._
(2018, 2020) and the results are reported in Table 2. We can see that by using
the PBE density and orbitals the trends do not change, i.e., the ordering of
the phases with (r)MGGAC is still correct, whereas it is not the case with
(r2)SCAN. There are some differences between the VASP self-consistent and
WIEN2k non-self-consistent results, however they are unimportant for the
conclusion. Overall, this indicates that the correct behavior of (r)MGGAC is
essentially functional-driven.

Next, we consider in detail the results for the structural, electronic, and
magnetic properties of $rs$-AF2 MnO, calculated with the semilocal and onsite
hybrid functionals. Our results are summarized in Table 3 and are compared to
experimental results (see Ref. Schrön  _et al._ , 2010 and references
therein). The experimental equilibrium lattice constant $a_{0}\sim 4.43$ Å is
best reproduced by PBE. (r2)SCAN and SCAN-$\alpha$ ($\alpha=0.25$) are also
pretty accurate. For the bulk modulus $B_{0}$, PBE is again the best method,
and PBE-$\alpha$ can be also very accurate depending on the value of $\alpha$,
while (r)MGGAC give values of $B_{0}$ that are too large by about 30 GPa.

The experimental band gap of 3.6$-$3.9 eV is strongly underestimated by at
least 2 eV by all methods in Table 3. In a recent study Tran  _et al._ (2018)
it was shown that among the fast semilocal DFT methods only the GLLB-SC Kuisma
_et al._ (2010) and Sloc Finzel and Baranov (2017) functionals are able to
give band gaps of MnO similar to experiment. Other studies Rödl  _et al._
(2009); Schiller  _et al._ (2015) have shown that $GW$ can be accurate
depending on the input orbitals, while DMC gives a band gap that is too large
by nearly 1 eV. The magnetic moment $\mu$ on the Mn atom was calculated inside
the atomic basin as defined by the quantum theory of atoms in molecules of
BaderBader (1990, 1991). The experimental value of 4.58 $\mu_{B}$ is best
reproduced by the MGGAC, rMGGAC, and PBE-$\alpha$($\alpha=0.25$) functionals
that give 4.52, 4.52, and 4.65 $\mu_{\text{B}}$, respectively.

Finally, we show in Fig. 4 the results for the other antiferromagnetic TMOs
considered in this work: FeO, CoO, and NiO. As for MnO, their ground-state
phase is $rs$-AF2 Roth (1958), and we can see that this is correctly predicted
by the MGGAC and rMGGAC functionals.

Neither experimental data nor values obtained from highly accurate methods
like DMC or RPA seem to be available for FeO. Therefore, a comparison of the
values of $\Delta E$ can be made only between DFT methods. As in the cases of
MnO and CoO (see below), (r)MGGAC lead to the largest positive values of
$\Delta E$, too. They are clearly larger than for all other methods. Next come
the onsite hybrid SCAN-$\alpha$ and $U$\- and/or vdW-corrected functionals
like SCAN+rVV10+$U$ or HSE+TS. Only PBE and SCAN lead erroneously to $zb$-AF1
as the ground-state phase, while a small amount of Hartree-Fock exchange
$\alpha$ or small value of $U$ is enough to get the correct ground-state phase
$rs$-AF2.

For CoO, the experimental and DMC values (see Ref. Saritas  _et al._ , 2018
and references therein) are around $400$ meV/f.u. and are reproduced very
accurately by (r)MGGAC, as visible on Fig. 4. Considering the onsite hybrid
functionals based on SCAN, using larger values of $\alpha$ leads to better
agreement with experiment and DMC, however the values are twice too small even
with $\alpha=0.55$. Values from the literature Schiller  _et al._ (2015); Peng
and Lany (2013); Peng and Perdew (2017); Saritas  _et al._ (2018) obtained
with other methods are also shown in Fig. 4. As for MnO, the popular methods
PBE, SCAN, and HSE06 fail since they predict $zb$-AF1 to be more stable than
$rs$-AF2, but adding a Hubbard $U$ and/or vdW correction helps to get the
correct trend.

As for FeO, no reference data is available for NiO. By inspecting the DFT
results, we can see that the situation is quite different compared to MnO,
FeO, and CoO. Firstly, all functionals, without exception, lead to $rs$-AF2 as
the ground-state phase. Secondly, the MGGAC and rMGGAC values are basically
the smallest in magnitude, while the reverse was obtained for MnO, FeO, and
CoO. The largest values of $\Delta E$ are obtained with the onsite hybrids
SCAN-$\alpha$ and PBE-$\alpha$, and are much larger than what is obtained with
the other methods. It should also be noted that the calculated energetic
differences between the $rs$-AF2 and $zb$-AF1 phases are much larger than for
the other studied TMOS by nearly one order of magnitude.

In order to show that the (r)MGGAC functionals lead to the correct energy
ordering also for other phases of the studied systems, we considered the
$wz$-AF1 phase Schrön  _et al._ (2010). The values of $\Delta E^{wz-
rs}=E_{wz}^{\text{AF1}}-E_{rs}^{\text{AF2}}$ are presented in Tables S2 and S3
of Ref. Ghosh  _et al._ , 2021b. For the four systems the (r)MGGAC values of
$\Delta E^{wz-rs}$ are positive, which is the correct trend. Reference
(experimental and DMC) values for $\Delta E^{wz-rs}$ are available only for
CoO Saritas  _et al._ (2018), and we can see that the magnitude of $\Delta
E^{wz-rs}$ obtained with (r)MGGAC is too large by roughly $\sim 200$ meV/f.u.
compared to the reference values.

At this point one may wonder why the (r)MGGAC functionals work much better
than PBE and (r2)SCAN for the relative stability of the phases of the studied
TMOs. In order to address this question, we performed calculations with the
WIEN2k code and decomposed $E_{\text{xc}}$ into contributions from the atomic
spheres and the interstitial. We show in Table 4 the total as well as the
decomposed values of $E_{\text{xc}}$ obtained from SCAN and MGGAC for the
$rs$-AF2 and $zb$-AF1 phases of MnO. MGGAC stabilizes $rs$-AF2 over $zb$-AF1
by 1515 meV/f.u., while SCAN gives a smaller value, 1324 meV/f.u., and this
difference of 191 meV/f.u. leads as discussed above to the correct ground
state with MGGAC but not with SCAN. By looking from which region this
difference comes, there is virtually no contribution from the O atom, a very
small one from Mn, while most of the stabilization comes from the interstitial
region. Of course the exact values depend weakly on the chosen atomic sphere
sizes (integration ranges). Therefore, we show in Fig. 5 the difference
$\epsilon_{\text{xc}}^{\text{SCAN}}-\epsilon_{\text{xc}}^{\text{MGGAC}}$
between the SCAN and MGGAC xc energy density [the integrand in Eq. (1)]. There
are large positive and negative differences inside the atomic sphere since the
corresponding density is large, but as shown in Table 4 these differences are
very similar for the two phases and cancel to a large extent. However, in the
interstitial regions we can clearly see that this difference is close to zero
in $zb$-AF1, but much more positive in $rs$-AF2. A positive difference
$\epsilon_{\text{xc}}^{\text{SCAN}}-\epsilon_{\text{xc}}^{\text{MGGAC}}$ means
a more negative integrand of $E_{\text{xc}}$ for MGGAC and thus a
stabilization of the $rs$-AF2 structure. An examination of the difference of
the enhancement factors
$F_{\text{xc}}^{\text{SCAN}}-F_{\text{xc}}^{\text{MGGAC}}$ (Fig. 6) shows
partially a similar picture (with opposite sign) as in Fig. 5, namely large
positive and negative differences inside the atomic spheres, which at the end
are similar for both $rs$-AF2 and $zb$-AF1 phases and thus cancel. On the
other hand this difference is more negative in the interstitial of the
$rs$-AF2 phase and together with a larger density gives the dominant
contribution.

Since now we know from which region (the interstitial) the effect comes, we
can mention what are the values of $r_{s}$, and the $\sigma$-spin $s_{\sigma}$
and $\alpha_{\sigma}^{\text{iso}}$ in this region (see Fig. 7). The value of
$r_{s}$ corresponding to the red interstitial region of the $rs$-AF2 phase in
$\epsilon_{\text{xc}}^{\text{SCAN}}-\epsilon_{\text{xc}}^{\text{MGGAC}}$ in
Fig. 5, or equivalently in blue in Fig. 6, is around 2 Bohr, while in the
large interstitial region of the $zb$-AF1 phase it is more around 3 Bohr,
which corresponds to lower density than in $rs$-AF2. $s_{\sigma}$ is quite
small (from 0 to 0.5), while $\alpha_{\sigma}^{\text{iso}}$ is moderately
large (from 2 to 3). However, it should be also noted that
$\alpha_{\sigma}^{\text{iso}}$ can reach values up to 8 in some regions just
outside of the $3d$ density (where also $s_{\sigma}$ increases up to 1.7) and
in the interstitial of $zb$-AF1 (where $s_{\sigma}$ is very small).

We can now compare the corresponding values of the enhancement factors
$F_{\text{xc}}$ of SCAN and MGGAC shown in Fig. 1. For the relevant values in
the interstitial region of $r_{s}$ (in the range 2-3), $s$ (below $\sim 0.7$),
and $\alpha_{\sigma}^{\text{iso}}$ (above $\sim 2$), we can see that
$F_{\text{xc}}^{\text{MGGAC}}$ is more positive than
$F_{\text{xc}}^{\text{SCAN}}$. Thus, since $\epsilon_{\text{xc}}$ is
proportional to $-n^{4/3}F_{\text{xc}}$ and the electron density $n$ in the
interstitial region is larger in $rs$-AF2 than in $zb$-AF1, then
$\epsilon_{\text{xc}}^{\text{MGGAC}}-\epsilon_{\text{xc}}^{\text{SCAN}}$ (and
$E_{\text{xc}}^{\text{MGGAC}}-E_{\text{xc}}^{\text{SCAN}}$, see Table 4) is
more negative in $rs$-AF2 than in $zb$-AF1, leading to the observed
stabilization of $rs$-AF2 with MGGAC.

##  IV Conclusions

In this work, the ability of various DFT methods to predict the correct
ground-state phase of antiferromagnetic transition-metal monoxides has been
studied. The case of MnO has been considered in more details. The conclusions
are the following. The meta-GGAs MGGAC and rMGGAC provide the correct energy
ordering of the $rs$-AF2 and $zb$-AF1 phases of MnO, while it is not the case
with the other popular functionals PBE, SCAN, and HSE06. Furthermore, the
relative energies of the two phases are in very nice agreement with the values
obtained with the highly accurate DMC and RPA methods. It should be underlined
that the very good (r)MGGAC results have been obtained without addition of vdW
or Hubbard $U$ correction. With other semilocal popular methods like PBE or
SCAN, it is necessary to add a $U$ or vdW correction to obtain the correct
ordering of the two phases. Since also the (r)MGGAC relative energy is very
accurate, this may indicate that (r)MGGAC is more accurate in the atomic
regions (where, alternatively, $U$ can also be added to improve the
description), but also in the interstitial region (where a vdW correction may
be helpful). We also showed in the present work that a relatively high amount
of Hartree-Fock exchange (which plays the same role as $U$) in onsite hybrids
is necessary to obtain the correct ordering of the phases of MnO.

Our results suggest that the semilocal (r)MGGAC functionals may be good
alternatives to the costly DMC and RPA methods to predict the ground-state
phase of strongly correlated systems. The low computational cost of semilocal
methods is certainly helpful for studying potentially interesting
functionalities of complex materials, such as heterostructure, cathode
materials, or alloys, where various structural phases may compete.

###  IV.1 Computational details

The calculations on MnO, FeO, CoO, and NiO with the r2SCAN and (r)MGGAC meta-
GGA functionals are performed using the Vienna Ab initio Simulation Package
(VASP) code Kresse and Hafner (1993); Kresse and Furthmüller (1996, 1996);
Kresse and Joubert (1999), which is based on the projector augmented wave
method Blöchl (1994); Kresse and Joubert (1999). The used pseudopotentials
correspond to valence electron configurations $3d^{6}4s^{1}$ for Mn,
$3d^{7}4s^{1}$ for Fe, $3d^{8}4s^{1}$ for Co, $3d^{9}4s^{1}$ for Ni, and
$2s^{2}p^{4}$ for O.

The calculations with the onsite hybrid functionals are done using the all-
electron WIEN2k code Blaha  _et al._ (2018, 2020), which is based on the
augmented plane wave plus local orbitals method Singh and Nordström (2006);
Karsai  _et al._ (2017).

All VASP and WIEN2k calculations are done with parameters (size of basis set,
k-mesh, etc.) that should lead to well converged results.

To calculate the equilibrium ground-state properties (total energy, volume,
isothermal bulk modulus and its pressure derivative) we fitted the total
energy versus volume data with the Birch-Murnaghan isothermal equation of
state.

###### Acknowledgements.

A.G. would like to thank INSPIRE fellowship, DST, India for financial support.
The first-principle calculations within VASP are performed using the high
performance computing (HPC) clusters of IITH, Hyderabad, KALINGA and NISERDFT
high performance computing (HPC) clusters of NISER, Bhubaneswar.

## References

  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Burke (2012) K. Burke, J. Chem. Phys. 136, 150901 (2012).
  * Jones (2015) R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
  * Becke (2014) A. D. Becke, J. Chem. Phys. 140, 18A301 (2014).
  * Perdew and Schmidt (2001) J. P. Perdew and K. Schmidt, AIP Conf. Proc. 577, 1 (2001).
  * Svane and Gunnarsson (1990) A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65, 1148 (1990).
  * Peng and Lany (2013) H. Peng and S. Lany, Phys. Rev. B 87, 174113 (2013).
  * Schrön  _et al._ (2010) A. Schrön, C. Rödl, and F. Bechstedt, Phys. Rev. B 82, 165109 (2010).
  * Schiller  _et al._ (2015) J. A. Schiller, L. K. Wagner, and E. Ertekin, Phys. Rev. B 92, 235209 (2015).
  * Peng and Perdew (2017) H. Peng and J. P. Perdew, Phys. Rev. B 96, 100101 (2017).
  * Kanan and Carter (2012) D. K. Kanan and E. A. Carter, J. Phys. Chem. C 116, 9876 (2012).
  * Toroker and Carter (2013) M. C. Toroker and E. A. Carter, J. Mater. Chem. A 1, 2474 (2013).
  * Peng and Lany (2012) H. Peng and S. Lany, Phys. Rev. B 85, 201202 (2012).
  * Gopal  _et al._ (2004) P. Gopal, N. A. Spaldin, and U. V. Waghmare, Phys. Rev. B 70, 205104 (2004).
  * Archer  _et al._ (2011) T. Archer, C. D. Pemmaraju, S. Sanvito, C. Franchini, J. He, A. Filippetti, P. Delugas, D. Puggioni, V. Fiorentini, R. Tiwari, and P. Majumdar, Phys. Rev. B 84, 115114 (2011).
  * Franchini  _et al._ (2005) C. Franchini, V. Bayer, R. Podloucky, J. Paier, and G. Kresse, Phys. Rev. B 72, 045132 (2005).
  * Roth (1958) W. L. Roth, Phys. Rev. 110, 1333 (1958).
  * Zhang  _et al._ (2020) Y. Zhang, J. Furness, R. Zhang, Z. Wang, A. Zunger, and J. Sun, Phys. Rev. B 102, 045112 (2020).
  * Sai Gautam and Carter (2018) G. Sai Gautam and E. A. Carter, Phys. Rev. Materials 2, 095401 (2018).
  * Tran  _et al._ (2018) F. Tran, S. Ehsan, and P. Blaha, Phys. Rev. Materials 2, 023802 (2018).
  * Tran  _et al._ (2020) F. Tran, G. Baudesson, J. Carrete, G. K. H. Madsen, P. Blaha, K. Schwarz, and D. J. Singh, Phys. Rev. B 102, 024407 (2020).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * Schrön  _et al._ (2012) A. Schrön, C. Rödl, and F. Bechstedt, Phys. Rev. B 86, 115134 (2012).
  * Tran  _et al._ (2006) F. Tran, P. Blaha, K. Schwarz, and P. Novák, Phys. Rev. B 74, 155108 (2006).
  * Marsman  _et al._ (2008) M. Marsman, J. Paier, A. Stroppa, and G. Kresse, J. Phys.: Condens. Matter 20, 064201 (2008).
  * Liu  _et al._ (2019) P. Liu, C. Franchini, M. Marsman, and G. Kresse, J. Phys.: Condens. Matter 32, 015502 (2019).
  * Jana  _et al._ (2018a) S. Jana, A. Patra, and P. Samal, J. Chem. Phys. 149, 094105 (2018a).
  * Jana and Samal (2019) S. Jana and P. Samal, Phys. Chem. Chem. Phys. 21, 3002 (2019).
  * Jana  _et al._ (2020a) S. Jana, A. Patra, L. A. Constantin, and P. Samal, J. Chem. Phys. 152, 044111 (2020a).
  * Jana  _et al._ (2020b) S. Jana, B. Patra, S. Śmiga, L. A. Constantin, and P. Samal, Phys. Rev. B 102, 155107 (2020b).
  * Solovyev and Terakura (1998) I. V. Solovyev and K. Terakura, Phys. Rev. B 58, 15496 (1998).
  * Engel and Schmid (2009) E. Engel and R. N. Schmid, Phys. Rev. Lett. 103, 036404 (2009).
  * Erten  _et al._ (2011) O. Erten, O. N. Meetei, A. Mukherjee, M. Randeria, N. Trivedi, and P. Woodward, Phys. Rev. Lett. 107, 257201 (2011).
  * Hossain  _et al._ (2021) M. Hossain, J. De, and J. Bhattacharjee, J. Phys. Chem. A 125, 6805 (2021).
  * Hossain and Bhattacharjee (2020) M. Hossain and J. Bhattacharjee, J. Chem. Phys. 153, 144103 (2020).
  * Massidda  _et al._ (1997) S. Massidda, A. Continenza, M. Posternak, and A. Baldereschi, Phys. Rev. B 55, 13494 (1997).
  * Kotani and van Schilfgaarde (2008) T. Kotani and M. van Schilfgaarde, J. Phys.: Condens. Matter 20, 295214 (2008).
  * Rödl  _et al._ (2009) C. Rödl, F. Fuchs, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 79, 235114 (2009).
  * Kuneš  _et al._ (2008) J. Kuneš, A. Lukoyanov, V. Anisimov, R. T. Scalettar, and W. E. Pickett, Nature Mater. 7, 198 (2008).
  * Mandal  _et al._ (2019) S. Mandal, K. Haule, K. Rabe, and K. M. Vanderbilt, Npj Comput. Mater. 5, 1 (2019).
  * Mukherjee  _et al._ (2014) A. Mukherjee, N. D. Patel, S. Dong, S. Johnston, A. Moreo, and E. Dagotto, Phys. Rev. B 90, 205133 (2014).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Sabatini  _et al._ (2013) R. Sabatini, T. Gorni, and S. de Gironcoli, Phys. Rev. B 87, 041108(R) (2013).
  * Jana  _et al._ (2018b) S. Jana, A. Patra, and P. Samal, J. Chem. Phys. 149, 044120 (2018b).
  * Patra  _et al._ (2020a) A. Patra, S. Jana, and P. Samal, J. Chem. Phys. 153, 184112 (2020a).
  * Jana  _et al._ (2021a) S. Jana, S. K. Behera, S. Śmiga, L. A. Constantin, and P. Samal, J. Chem. Phys. 155, 024103 (2021a).
  * Jana  _et al._ (2020c) S. Jana, L. A. Constantin, and P. Samal, J. Chem. Theory Comput. 16, 974 (2020c).
  * Jana  _et al._ (2020d) S. Jana, A. Patra, S. Śmiga, L. A. Constantin, and P. Samal, J. Chem. Phys. 153, 214116 (2020d).
  * Ghosh  _et al._ (2021a) A. Ghosh, S. Jana, M. Niranjan, S. K. Behera, L. A. Constantin, and P. Samal, Journal of Physics: Condensed Matter (2021a).
  * Tran  _et al._ (2016) F. Tran, J. Stelzl, and P. Blaha, J. Chem. Phys. 144, 204120 (2016).
  * Sengupta  _et al._ (2018) N. Sengupta, J. E. Bates, and A. Ruzsinszky, Phys. Rev. B 97, 235136 (2018).
  * Shahi  _et al._ (2018) C. Shahi, J. Sun, and J. P. Perdew, Phys. Rev. B 97, 094111 (2018).
  * Peng  _et al._ (2016) H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Phys. Rev. X 6, 041005 (2016).
  * Mo  _et al._ (2017) Y. Mo, R. Car, V. N. Staroverov, G. E. Scuseria, and J. Tao, Phys. Rev. B 95, 035118 (2017).
  * Jana  _et al._ (2018c) S. Jana, K. Sharma, and P. Samal, J. Chem. Phys. 149, 164703 (2018c).
  * Jana  _et al._ (2019) S. Jana, K. Sharma, and P. Samal, J. Phys. Chem. A 123, 6356 (2019).
  * Patra  _et al._ (2019a) B. Patra, S. Jana, L. A. Constantin, and P. Samal, Phys. Rev. B 100, 155140 (2019a).
  * Patra  _et al._ (2019b) A. Patra, S. Jana, and P. Samal, J. Phys. Chem. A 123, 10582 (2019b).
  * Patra  _et al._ (2021a) B. Patra, S. Jana, L. A. Constantin, and P. Samal, J. Phys. Chem. C 125, 4284 (2021a).
  * Jana  _et al._ (2021b) S. Jana, S. K. Behera, S. Śmiga, L. A. Constantin, and P. Samal, New J. Phys. (2021b).
  * Patra  _et al._ (2020b) A. Patra, B. Patra, L. A. Constantin, and P. Samal, Phys. Rev. B 102, 045135 (2020b).
  * Patra  _et al._ (2021b) A. Patra, S. Jana, P. Samal, F. Tran, L. Kalantari, J. Doumont, and P. Blaha, J. Phys. Chem. C 125, 11206 (2021b).
  * Tran  _et al._ (2021) F. Tran, J. Doumont, L. Kalantari, P. Blaha, T. Rauch, P. Borlido, S. Botti, M. A. L. Marques, A. Patra, S. Jana, and P. Samal, J. Chem. Phys. 155, 104103 (2021), https://doi.org/10.1063/5.0059036 .
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and Y. Wang, Phys. Rev. B 54, 16533 (1996).
  * Furness  _et al._ (2020) J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020).
  * Seidl  _et al._ (1996) A. Seidl, A. Görling, P. Vogl, J. A. Majewski, and M. Levy, Phys. Rev. B 53, 3764 (1996).
  * Yang  _et al._ (2016) Z.-h. Yang, H. Peng, J. Sun, and J. P. Perdew, Phys. Rev. B 93, 205205 (2016).
  * Perdew  _et al._ (2017) J. P. Perdew, W. Yang, K. Burke, Z. Yang, E. K. U. Gross, M. Scheffler, G. E. Scuseria, T. M. Henderson, I. Y. Zhang, A. Ruzsinszky, H. Peng, J. Sun, E. Trushin, and A. Görling, Proc. Natl. Acad. Sci. U. S. A. 114, 2801 (2017).
  * Sun  _et al._ (2013) J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P. Perdew, Phys. Rev. Lett. 111, 106401 (2013).
  * Della Sala  _et al._ (2016) F. Della Sala, E. Fabiano, and L. A. Constantin, Int. J. Quantum Chem. 22, 1641 (2016).
  * Becke and Roussel (1989) A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761 (1989).
  * Patra  _et al._ (2019c) B. Patra, S. Jana, L. A. Constantin, and P. Samal, Phys. Rev. B 100, 045147 (2019c).
  * Perdew  _et al._ (2014) J. P. Perdew, A. Ruzsinszky, J. Sun, and K. Burke, J. Chem. Phys. 140, 18A533 (2014).
  * Aschebrock and Kümmel (2019) T. Aschebrock and S. Kümmel, Phys. Rev. Res. 1, 033082 (2019).
  * Novák  _et al._ (2006) P. Novák, J. Kuneš, L. Chaput, and W. E. Pickett, Phys. Status Solidi B 243, 563 (2006).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  * Vydrov and Van Voorhis (2010) O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 133, 244103 (2010).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Kim  _et al._ (2013) M.-C. Kim, E. Sim, and K. Burke, Phys. Rev. Lett. 111, 073003 (2013).
  * Blaha  _et al._ (2018) P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, J. Luitz, R. Laskowski, F. Tran, and L. D. Marks,  _WIEN2k: An Augmented Plane Wave plus Local Orbitals Program for Calculating Crystal Properties_ (Vienna University of Technology, Austria, 2018).
  * Blaha  _et al._ (2020) P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, J. Chem. Phys. 152, 074101 (2020).
  * Kuisma  _et al._ (2010) M. Kuisma, J. Ojanen, J. Enkovaara, and T. T. Rantala, Phys. Rev. B 82, 115106 (2010).
  * Finzel and Baranov (2017) K. Finzel and A. I. Baranov, Int. J. Quantum Chem. 117, 40 (2017).
  * Bader (1990) R. F. W. Bader,  _Atoms in Molecules: A Quantum Theory_ (Oxford University Press, Oxford, 1990).
  * Bader (1991) R. F. W. Bader, Chem. Rev. 91, 893 (1991).
  * Saritas  _et al._ (2018) K. Saritas, J. T. Krogel, and F. A. Reboredo, Phys. Rev. B 98, 155130 (2018).
  * Ghosh  _et al._ (2021b) A. Ghosh, S. Jana, M. K. Niranjan, T. Fabien, D. Wimberger, B. Peter, L. A. Constantin, and P. Samal, supporting information (2021b).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Singh and Nordström (2006) D. J. Singh and L. Nordström,  _Planewaves, Pseudopotentials, and the LAPW Method, 2nd ed._ (Springer, New York, 2006).
  * Karsai  _et al._ (2017) F. Karsai, F. Tran, and P. Blaha, Comput. Phys. Commun. 220, 230 (2017).


**footnotetext: Equal contribution11institutetext: Skoltech, Moscow, Russia
22institutetext: AIRI, Moscow, Russia 33institutetext: University of Hamburg,
Hamburg, Germany

# Label Attention Network for sequential multi-label classification: you were
looking at a wrong self-attention

Elizaveta Kovtun∗ 11 Galina Boeva∗ 11 Artem Zabolotnyi 11 Evgeny Burnaev 11 2
2 Martin Spindler 33 Alexey Zaytsev 11

###### Abstract

Most of the available user information can be represented as a sequence of
timestamped events. Each event is assigned a set of categorical labels whose
future structure is of great interest. For instance, our goal is to predict a
group of items in the next customer’s purchase or tomorrow’s client
transactions. This is a multi-label classification problem for sequential
data. Modern approaches focus on transformer architecture for sequential data
introducing self-attention for the elements in a sequence. In that case, we
take into account events’ time interactions but lose information on label
inter-dependencies. Motivated by this shortcoming, we propose leveraging a
self-attention mechanism over labels preceding the predicted step. As our
approach is a Label-Attention NETwork, we call it LANET. Experimental evidence
suggests that LANET outperforms the established models’ performance and
greatly captures interconnections between labels. For example, the micro-AUC
of our approach is 0.9536 compared to 0.7501 for a vanilla transformer. We
provide an implementation of LANET to facilitate its wider usage.

###### Keywords:

Multi-label classification, attention mechanism, sequential data

##  1 Introduction

The multi-label classification is a more natural setting than a binary or
multiclass classification since everything that surrounds us in the real world
is usually described with multiple labels [19]. The same logic can be
transferred to the sequence of timestamped events. Events in a sequence are
prone to be characterized by several categorical values instead of one mark.
There are numerous approaches to deal with the multi-label classification in
computer vision [6], natural language processing [37], or classic tabular data
framework [29]. However, the multi-label problem statement for event sequences
tends to receive less attention. So, we aim to confront such a lack of focus
and solve the labels’ set prediction problem for sequential timestamped data.
In short, it is a multi-label classification for sequential data.

Importantly, the model should predict a set of labels that correspond to the
next step by taking into account the content of previous label groups for a
sequence of events related to an object. The problem at hand is illustrated by
an example in Figure 1.

ID | Date | Labels | Features
---|---|---|---
1 | 12-08-2016 | [1, 18, 89] | [2.1, 0.4, 0.7]
| 21-08-2016 | [3, 8] | [0.3, 1.5]
| 23-08-2016 | [1, 18] | [1.7, 0.5]
| 28-08-2016 | [?, ?, …] |

Figure 1: A dataset sample of timestamped events that can be used within the
multi-label classification task and its visualization. ID refers to a
considered client identifier, while features are additional characteristics of
labels in a set. For example, if a set of labels are a group of transactions,
then features represent the amounts of particular transactions. Our model
should predict labels for the $28$th of August, given the history of previous
label sets. There are multiple labels to predict, so it is a multi-label
classification problem. The right figure provides a visual representation of
the problem in the left table.

The interaction between an object’s states at different timestamps is
essential for solving tasks with sequential data [10]. Therefore, expressive
and powerful models should be able to learn such interactions. Several neural
network architectures, such as transformers or recurrent neural networks, are
capable of doing this. For example, a transformer directly defines an
attention mechanism that measures how different timestamps in a sequence are
connected. However, the applications of modern deep learning methods are
limited [43], and they primarily focus on predicting labels for a sequence in
general. We refer to the graph of connections between states of an object at
different timestamps as a timestamp interaction graph.

Another connection worth exploring is the connection between different labels
and a need to consider the correlation between them [9]. This capability is
not inherent for the majority of models. With a large number of possible
labels, this problem becomes close to a recommender system problem. We name
the graph of connections between different labels a label interaction graph.

For single label-event data by using a position encoder, we take into account
both interaction between labels and timestamped events [16, 40]. For
sequential multi-label prediction, simultaneous consideration of both
timestamp interaction graph and label interaction graph can be crucial.
Typically, articles explore a particular side of the dependence that can be
explainable by domain bias. In sequential recommender systems, there is a
focus on connections between labels [23] with the incorporation of
convolutional neural networks [28] as well as the attention mechanism [42].
Direct models for event sequences [11] also prefer identification of
interactions between timestamps [44]. So, in time series and natural language
processing, it is more common to handle a timestamp interaction graph via an
attention mechanism or recurrent neural networks.

|  LSTM |  TransformerBase |  CLASS2C2AE |  LANET (ours)
---|---|---|---|---
Micro-AUC |  2 |  3.2 |  3.6 |  1.2
Macro-AUC |  2 |  3.2 |  3.6 |  1.2
Micro-F1 |  1.8 |  3.4 |  3.4 |  1.6
Macro-F1 |  1.8 |  3.2 |  2.8 |  2.2
Table 1: Mean rank for different metrics averaged over $5$ considered
datasets. We want to minimize rank, as the best method has rank $1$. The best
values are in bold and the second best values are underlined.

Our LANET aims at recovering label interaction graphs, as we believe it is
more vital in many problems. The algorithm aggregates past information in a
specially constructed set of label representations. The set serves as input to
a transformer encoder. The transformer encoder updates label embeddings via
self-attention, uncovering their interactions. Finally, we predict a set of
labels for the current model based on the model’s output that encompasses
knowledge of label relationships. Moreover, we can process long sequences this
way, as now the attention evaluation is quadratic in the number of labels, not
the sequence length. As our model does not use external information, it
produces self-supervised embeddings.

#### Contributions.

We have developed a transformer-based architecture with self-attention between
labels to deal with multi-label classification for event sequences. Our main
contributions are the following:

  * •

We introduce LANET architecture for the prediction of a label set for the
current event taking the information from previous events. The architecture’s
peculiarity is a calculation of self-attention between label views. The scheme
of our approach is presented in Figure 2.

  * •

LANET outperforms transformer-based models that focus on self-attention
between timestamps. We evaluate all metrics on various datasets to prove the
insufficiency of a common sequential placement-oriented perspective. See Table
1 for a high-level comparison of different approaches. Furthermore, we provide
an ablation study of the proposed architecture.

  * •

Due to LANET structure, it efficiently captures label interactions via the
attention matrix.

Figure 2: LANET architecture for multi-label classification of event sequences
with a total number of labels $K$. The model input comprises the concatenation
of the object’s label, time, feature, and ID embeddings. After self-attention
layers, the model outputs a vector of confidence scores, where each element is
associated with the probability of the presence of a particular label in the
next event. The self-attention layer in LANET acts between labels, not
temporal representations.

##  2 Related Work

The multi-label classification problem statement emerges in many diverse
domains, e.g., text categorization or image tagging, all of which entail their
own peculiarities and challenges. The review [41] explores foundations in
multi-label learning, discussing the well-established methods as well as the
most recent approaches. Emerging trends are covered in a more recent review
[19].

Firstly, we examine loss functions tailored for the multi-label setting and
into some methods for composing label set prediction. Secondly, we overview
the usage of RNNs in the multi-label classification task. Thirdly, we review
how to capture the label dependencies. Then, we discuss an association with a
sequential recommendation problem. Lastly, we summarise the literature
analysis and present the identified research gap.

#### Loss functions and ways for label set composition in multi-label problem.

The paper [21] studies the theoretical background for main approaches to the
reduction of multi-label classification problem to a series of binary or
multi-class problems. In particular, they show that considered reductions
implicitly optimize for either Precision@k or Recall@k. The choice of the
correct reduction should be based on the ultimate performance measure of
interest. In [17], the authors propose an improved loss function for pairwise
ranking in multi-label image classification task that is easier to optimize.
Also, they discuss an approach based on the estimation of the optimal
confidence thresholds for the label decision part of the model that determines
which labels to include in the final prediction. The task of multi-label text
classification is the topic of [5]. The authors construct an end-to-end deep
learning framework called ML-Net. ML-Net consists of two parts: a label
prediction network and a label count prediction network. In order to get the
final set of labels, confidence scores generated from the label prediction
network are ranked, and then the top $K_{top}$ labels are predicted. A
separate label count network predicts $K_{top}$.

#### Neural networks for multi-label classification.

In [38], the authors use the RNN model to solve a multi-label classification
problem. They address the issue that RNNs produce sequential outputs, so the
target labels should be ordered. The authors propose to dynamically order the
ground truth labels based on the model predictions, which contributes to
faster training and alleviates the effect of duplicate generation. In turn,
[31] consider the transformation of a multi-label classification problem to a
sequence prediction problem with an RNN decoder. They propose a new learning
algorithm for RNN-based decoders which does not rely on a predefined label
order. Consequently, the model explores diverse label combinations,
alleviating the exposure bias. The work [26] examines the same problem
statement of multi-label classification in an event stream as we do. The
authors’ model targets capturing temporal and probabilistic dependencies
between concurrent event types by encoding historical information with a
transformer and then leveraging a conditional mixture of Bernoulli experts.

#### Approaches to leveraging label dependencies.

The authors of [15] construct a model called C-Tran for a multi-label image
classification task that leverages Transformer architecture that encourages
capturing the dependencies among image features and target labels. The key
idea is to train the model with label masking. The authors in [39] propose DNN
architecture for solving the multi-label classification task, which
incorporates the construction of label embeddings with feature and label
interdependency awareness. A label-correlation sensitive loss improves the
efficiency of the constructed model. Another popular way to consider label
relationships is to use Graph Neural Networks as a part of the pipeline.
Namely, [22] captures the correlation between the labels in the task of Multi-
Label Text Classification by adopting Graph Attention Network (GAT). They
predict the final set of labels combining feature vectors from BiLSTM and
attended label features from GAT. Event sequences processing also tries to
derive dependencies between different event types and consider specific
attention mechanisms [20]. The closest to our LANET approach [18] in the
literature explores the relationship between different series for multivariate
time series classification. The authors propose using attention in step-wise
and channel-wise fashion to produce embeddings, which are then forwarded by a
classification head.

#### Sequential recommender systems.

Another close neighbour of our problem statement is the problem of sequential
recommender system construction [34, 23]. In this case, we have a lot of
possible labels and should sort them with respect to the probability of
occurrence next time. Typically, the estimation of embeddings for all possible
labels/items is a part of a pipeline. Existing approaches use neural networks
for sequential data such as LSTM [35] as well as attention mechanism [33]. We
want to highlight the statement related to the usage of only recent past data
for prediction [16]. However, millions of possible labels typically lead to
more classic techniques in this area with specific loss functions and methods.

#### Research gap.

Almost all existing methods are applied to multi-label learning for images,
natural language texts or time series. However, there is a need for approaches
to the multi-label classification of event sequences. Such characteristics of
events as time-dependency and heterogeneous attributes introduce additional
difficulties as well as opportunities for their processing and learning of
label correlations. Notably, such data structure encourages the usage of the
attention mechanism, as the sets of possible labels in the past and the future
coincide, and it is natural to leverage such a connection in an architecture.

##  3 Methods

###  3.1 Multi-label sequential classification

We consider a multi-label classification for a sequence
$S=\\{(X_{i},Y_{i})\\}_{i=1}^{t-1}$. It consists of a set of labels $Y_{i}$
and a set of features $X_{i}$ specific for each timestamp from $1$ to $t-1$.
The index corresponds to the time of the event, so $(X_{1},Y_{1})$ is the
information about the first event, and $(X_{t-1},Y_{t-1})$ is the information
about the last observed event. A set $Y_{i}\subseteq\mathcal{Y}$, where
$\mathcal{Y}=\\{1,2,\dots,K\\}$ is the set of all possible labels. A set size
of $X_{i}$ is equal to a set size of $Y_{i}$. Each label from $Y_{i}$ is
entailed with the numerical or categorical feature from $X_{i}$ in the
corresponding position. We also can have an additional feature vector
$\mathbf{z}$ describing a considered sequence $S$ as a whole, e.g. user ID.
The goal of sequential multi-label classification is to predict a set of
labels $Y_{t}$ for the next timestamp.

We construct a function $f(\cdot)\in[0,1]^{K}$ that takes the historical
information on events as input and outputs the vector of scores for each of
$K$ labels. Each score is an estimate of the probability of a label to be
present in the next event-related label set. In our setting, we limit the size
of the past available to our model.
$S^{t}=\\{(X_{j},Y_{j})\\}_{j=t-\tau}^{t-1}$, where $\tau$ means a number of
events preceding the considered event with the timestamp $t$, which is
attributed with a target label set $Y_{t}$. So, more formally $f(\cdot)$ has
the form:

$f(X_{t-\tau},\ldots,X_{t-1},Y_{t-\tau},\ldots,Y_{t-1},\textbf{z})\in[0,1]^{K}$

to prepare scores for the prediction of $Y_{t}$.

To complete a prediction, we need a separate label decision model
$g(f(\cdot))$, that transforms confidence scores into labels. For example, we
compare the score for $k$-th label to a selected threshold $\beta_{k}$: if
$f_{k}(\cdot)>\beta_{k},k=1,\dots,K$, then the model predicts that the $k$-th
label is present. Thus, the model $g$ produces a label set
$\hat{Y}_{t}\subseteq\mathcal{Y}$ from the input confidence scores.

The resulting quality of the model depends both on the method $g(\cdot)$ for
selecting labels to a final set and the performance of $f(\cdot)$ to produce
confidence scores while we focus on working with $f(\cdot)$.

###  3.2 Dataset structure

We consider problems with the following dataset structure. See Figure 1 for an
example of a data sample.

  * •

ID. Each sequence that corresponds to a client or an organization has a unique
identifier. So in our case, the feature vector for a sequence description is
just a scalar value $z$. There may be thousands of unique IDs in the dataset.
Nevertheless, we as well provide the results of experiments without the usage
of ID information.

  * •

Date. Each event is associated with a date.

  * •

Labels. The label’s column provides a set of categorical variables $Y_{i}$
that describes an event. It might be a sequence of codes of transactions made
by a client in a particular day or the categories of the sold items.

  * •

Features. Each event is also characterized by a feature set $X_{i}$. It is
represented as a set of categorical or numerical values. Each component in a
set describes the corresponding label from a label set. Values of a feature
set may have a meaning of the sums of paid money per transaction code or the
number of the sold items per category.

###  3.3 Sample structure

We work within a sequential multi-label classification problem in the panel
setting. The considered datasets have the structure
$\mathcal{D}=\\{(S_{j},z_{j})\\}_{j=1}^{n}$. Each $S_{j}$ is a sequence of
pairs of feature sets and corresponding multi-label sets. We need to split
sequences $S_{j},j=1,\dots,n$ into samples for model training and testing. The
part of a sample that is passed to the model’s input $f(\cdot)$ incorporates
information on $\tau$ consecutive events from $S_{j}$ preceding the considered
timestamp $t$. The sample’s target is a label set of the event from $S_{j}$
with the considered timestamp $t$. Thus, to construct train and test samples
from the initial sequence $S_{j}$, we need to apply a sliding window approach.
Importantly, train and test datasets that share the samples related to the
same $z_{j},j=1,\dots,n$ organized in a non-overlapping and chronological way.

###  3.4 Our LANET approach

Most of the transformer-related models used for sequential multi-label
prediction use self-attention computation between consecutive input timestamps
representations. The LANET instead uses the self-attention between label
representations. So, it has the input that consists of $K$ vectors. Below, we
describe how to aggregate a sequence of size $\tau$ to $K$ vectors via an
Embedding layer. Then we define the Self-attention layer. To get the
predictions, we apply a Prediction layer.

#### Embedding layer.

We use the following approach to use different parts of input data for multi-
label event sequences:

  * •

_ID embedding:_ For IDs we learn an embedding matrix;

  * •

_Time embedding:_ For each timestamp, we know the value of $dt$, which is
equal to the difference in days between the considered and the previous
timestamp. We train an embedding for each observed $dt$ value. We also take
into account the order of events, so we look at the embeddings for positions:
$1,\dots,$ $\tau$ to add them to the $dt$ representation of the corresponding
timestamp;

  * •

_Amount embedding:_ We transform all amounts into bins, breaking down the
continuous sums of amounts on the interval. Each interval is assigned a unique
number. Then, for each unique number, we create embedding. A vector of
representations is learned similar to [8].

  * •

_Transformer label encoder:_ To construct the LANET’s input, we use the data
associated with a specific ID for $\tau$ timestamps. We concatenate the
representation of the labels that are encountered during $\tau$ timestamps
with the corresponding time and amount embeddings. If a label does not belong
to a history of $\tau$ last timestamps, we add vectors of zeros to it as time
and amount views. If a particular label occurs several times during previous
$\tau$ steps, then we construct time embeddings for each individual occurrence
and then sum them up to get the ultimate time view for that label. We do the
same when constructing amount embedding in the situation of label
reoccurrence.

So, as a result of the embedding layer, we have $K+1$ embedding vectors. The
first vector $\mathbf{e}_{0}$ corresponds to the embeddings of general (ID)
features for a sequence. All other vectors $\mathbf{e}_{k},k=1,\dots,K$ are
concatenation of the label embedding with the corresponding time and amount
embeddings. It turns out that vectors for not historically involved labels are
just label embeddings, as we are summing them up with zero vectors of time and
and amount views. Under the training of all embedding weights, they are
initialized from normal distribution $\mathcal{N}(0,1)$ and then optimized.

#### Self-attention layer.

After getting representations from our data, we move on to the composition of
our architecture. Let $E=\\{\mathbf{e}_{0},\dots,\mathbf{e}_{K}\\}$,
$\mathbf{e}_{i}\in\mathbb{R}^{d}$, be the sequence of input embeddings to the
Transformer encoder, where $\mathbf{e}_{0}$ corresponds to the embedding of
$\mathbf{z}$ and all other $\mathbf{e}_{i}$ correspond to embeddings captured
historical information from a label point of view. In Transformer
architecture, the influence of embedding $\mathbf{e}_{j}$ on embedding
$\mathbf{e}_{i}$ is obtained via self-attention. The attention weight
$\alpha_{ij}$ and an embedding update $\mathbf{e}_{i}^{\prime}$ are calculated
as:

$\mathbf{e}_{i}^{\prime}=\sum_{j=0}^{K}\alpha_{ij}(W^{v}\mathbf{e}_{j});\hskip
8.53581pt\alpha_{ij}=\mathrm{softmax}\left((W^{q}\mathbf{e}_{i})^{T}(W^{k}\mathbf{e}_{j})/\sqrt{d}\right),$

where $W^{k}$ is the key weight matrix, $W^{q}$ is the query weight matrix,
and $W^{v}$ is the value weight matrix. Such a procedure of embedding updates
can be repeated several times by specifying the number of Transformer encoder
layers.

#### Prediction layer.

The updated embeddings are processed with one fully-connected linear layer to
get final embeddings for each label
$\\{\mathbf{e}^{(final)}_{j}\\}_{j=1}^{K}$. In this way, we obtain
$\\{f_{j}\\}_{j=1}^{K}$ that are used in a multi-label classifier with
threshold $t_{j}$ selected separately for each label using a validation
sample.

##  4 Experiments

In this section, we present comparison of our approach with others and
relevant ablation studies. The code is available at GitHub
repository***https://github.com/E-Kovtun/order_prediction/tree/common_model.

###  4.1 Datasets

|  Dataset |  # events |  Median |  Max |  # unique |  Diff
---|---|---|---|---|---|---
|  |  |  set size |  set size |  labels |
|  Sales |  47 217 |  16 |  48 |  84 |  0.0632
|  Demand |  5 912 |  13 |  24 |  33 |  0.0957
|  Liquor |  291 029 |  14 |  66 |  107 |  0.0413
|  Transactions |  784 520 |  3 |  23 |  77 |  0.1079
|  Orders |  226 522 |  2 |  13 |  61 |  0.0518
Table 2: Characteristics of the datasets used in sequential multi-label
classification problems.

For a more versatile study of the algorithm we built, we conducted experiments
on five diverse datasets. Sales dataset [4] is historical sales data from
different shops. The labels relate to the item categories, and the amount is
the number of sold items for a particular category. Demand dataset [7]
describes historical product demand of several warehouses. The label feature
means the product category, and the amount feature relates to the
corresponding demand. Liquor dataset [27] presents information on the liquor
sales of the different stores. Each liquor label is identified by its
category. The amount of sold liquor is measured in litres. Transactions
dataset [2] contains histories of clients transactions. Each transaction is
described with the special MCC code and an involved amount of money. In the
original paper, the dataset has the name Gender. To demonstrate the
applicability of the proposed approach in the industry, we also consider a
private Orders dataset. It provides information on the restaurant orders of
different materials for beer production. The amount feature is related to the
volume of ordered materials in hectolitres.

The overall statistics for the considered datasets are given in Table 2. We
present the number of observed events, Median set size of all available label
sets $\mathrm{median}({|Y_{ij}|}_{i,j=1,1}^{n,t_{i}})$, Maximum size of the
label set that is encountered in the dataset
$\max({|Y_{ij}|}_{i,j=1,1}^{n,t_{i}})$, the number of unique labels $K$, and
Diff. Diff measures the imbalance of labels: we calculate 5$\%$ and 95$\%$
quantiles for frequencies of labels and take the difference between them. The
representation of labels is imbalanced in the majority of datasets, but we use
metrics that take into account this effect. All datasets are transformed to a
general structure described in Figure 1.

Dataset | Model |  Micro-AUC$\uparrow$ |  Macro-AUC$\uparrow$ |  Micro-F1$\uparrow$ |  Macro-F1$\uparrow$ |
---|---|---|---|---|---|---
| LSTM | 0.8670 | 0.7346 | 0.5600 | 0.4389 |
| TransformerBase | 0.8604 | 0.7206 | 0.5348 | 0.3751 |
| CLASS2C2AE | 0.8528 | 0.6881 | 0.5116 | 0.4217 |
Sales | LANET (ours) | 0.9069 | 0.7627 | 0.6235 | 0.4901 |
| LSTM | 0.8829 | 0.7633 | 0.6746 | 0.5929 |
| TransformerBase | 0.8624 | 0.7240 | 0.6491 | 0.5678 |
| CLASS2C2AE | 0.8342 | 0.7079 | 0.6738 | 0.5581 |
Demand | LANET (ours) | 0.8806 | 0.7373 | 0.7038 | 0.5908 |
| LSTM | 0.7520 | 0.7253 | 0.4729 | 0.2218 |
| TransformerBase | 0.7501 | 0.7197 | 0.4674 | 0.2279 |
| CLASS2C2AE | 0.9442 | 0.8058 | 0.4088 | 0.3638 |
Liquor | LANET (ours) | 0.9536 | 0.8728 | 0.7089 | 0.4508 |
| LSTM | 0.9719 | 0.8189 | 0.2133 | 0.2349 |
| TransformerBase | 0.9714 | 0.8098 | 0.1982 | 0.2285 |
| CLASS2C2AE | 0.9347 | 0.7324 | 0.1131 | 0.1501 |
Transactions | LANET (ours) | 0.9743 | 0.8328 | 0.2373 | 0.1067 |
| LSTM | 0.9800 | 0.9562 | 0.6827 | 0.6280 |
| TransformerBase | 0.9601 | 0.9316 | 0.6302 | 0.5631 |
| CLASS2C2AE | 0.9218 | 0.8964 | 0.6324 | 0.5979 |
Orders | LANET (ours) | 0.9847 | 0.9665 | 0.5773 | 0.5856 |
Table 3: Comparison of our method LANET with the baselines on five different
datasets. Best values are highlighted, and second-best values are underlined.

###  4.2 Baselines and training details

TransformerBase utilizes the transformer architecture [32]. Similar to many
previous works, it considers the self-attention between input vectors that
correspond to consecutive points in a sequence. We also note that it can be
viewed as a simplification of next-item-recommendation models [16]. _The
amount vector for each timestamp_ is the input vector, whose dimension is
equal to the number of the unique label values. Each position of the vector
contains the corresponding amount value. There is no need for label embedding,
as the information on historically involved labels is introduced through the
amount vector. Time and ID embeddings are constructed in the same way as in
LANET. We also use two self-attention layers. The output of the transformer
encoder is an updated time representation. Subsequently, Mean Pooling is
applied to this output to obtain the final vector of confidence estimates,
based on which we make a prediction on a set of labels for the current
timestamp.

We also use two RNN-based approaches. RNN (LSTM) replaces the transformer
encoder with LSTM-type encoder [13]. The model takes one by one point in time
and updates the hidden representation, which contains information about
previous elements. CLASS2C2AE is a model for multilabel-classification
presented in [39]. This model combines LSTM and autoencoder to improve its
performance.

For all models, we use the same data input given in Subsection 3.2. The
following manually selected hyperparameters were used for all experiments:
dimension of the embedding space is 128, batch size is 32. We use Adam
optimizer with learning rate initialized to $0.001$ and reduced via reduce-LR-
on-plateau scheduler.

#### LANET training details.

Our LANET model consists of two Transformer encoder layers with four heads in
each layer. Self-attention and layer implementation are inherited from the
basic PyTorch Transformer layer [1]. After receiving the updated vector
representations, we run them through the dropout layer, which is initialized
with the probability $0.3$. In the end, a fully-connected linear layer is
applied to the obtained vectors to get final label scores. LANET is trained
with the cross-entropy loss adapted for the multi-label task through
independent consideration of each label score.

#### Validation procedure.

The original datasets are divided into train, validation, and test sets. For
each ID, some part of the consecutive events go to the train set, and the rest
is allocated for validation and test. Importantly, the time periods in the
train, valid, and test sets do not overlap and are arranged in the order of
increasing time. We take 70% of the data samples for model training, 17% for
validation, and 13% for testing. All experiments were also carried out for
five random seeds, and the results after the experiments were averaged.

The article [36] provides an overview of metrics for multi-label
classification. We consider the most frequently used and comprehensive metrics
for the estimation of multi-label classification quality: Micro-AUC, Macro-
AUC, Micro-F1, and Macro-F1.

###  4.3 Main results

For all datasets, the models are trained with five different random seeds to
understand the value of results deviation from the result depending on the
selected random state. The standard deviation influences the metrics in the
third decimal place, while the superiority in performance is observed in the
first or second decimal place. To make our tables more readable, we omit the
standard deviation values. To calculate micro-F1 and macro-F1, we search for
thresholds to compare them with the obtained confidence scores and compose the
final label set. Thresholds are selected on the validation set by optimizing
F1-score for each label separately.

Table 3 presents the sequential multi-label classification metrics for
considered approaches. Our LANET gives the best results for diverse datasets.
The LSTM baseline is also suitable for capturing the label relationships in
general. We denote that it gives better results than the used version of
transformer. The reason can be the limited amount of data for training
available, so there are too many parameters for a transformer with timestamp-
based self-attention. Besides, we test the model suggested in [26] for our
problems. However, the results were not satisfying. We associate it with a
more complex structure and, consequently, greater sensitivity to model
parameters.

###  4.4 Ablation study

In this section, we provide results of ablation studies. The results are for
Demand dataset, if not mentioned otherwise, as for most of the other datasets
conclusions are similar.

#### Usage of both label-attention and timestamp-attention.

It is worth exploring what perspective of self-attention calculation is more
critical: between labels or between timestamps. Table 4 answers this question.
Label-attention is our basic implementation. Time-attention is a case when we
consider attention only between timestamp views. Concat-attention implies
getting confidence scores from concatenation of label-attention and time-
attention. We learn the importance coefficients of two forms of attentions and
use them as weights when summing attention views in the case of gated-
attention. Absence-indication experiment consists in adding a learnable vector
to the input embeddings if the particular label is not involved in the
considered history. As a result, Label-attention is the most beneficial, and
the inclusion of both types of attention decreases the quality.

|  Micro-AUC$\uparrow$ |  Macro-AUC$\uparrow$ |  Micro-F1$\uparrow$ |  Macro-F1$\uparrow$
---|---|---|---|---
Label-attention |  0.881 $\pm$ 0.007 |  0.737 $\pm$ 0.017 |  0.704 $\pm$ 0.018 |  0.591 $\pm$ 0.003
Time-attention |  0.837 $\pm$ 0.004 |  0.682 $\pm$ 0.007 |  0.674 $\pm$ 0.009 |  0.588 $\pm$ 0.002
Concat-attention |  0.835 $\pm$ 0.001 |  0.681 $\pm$ 0.010 |  0.666 $\pm$ 0.000 |  0.587 $\pm$ 0.001
Gated-attention |  0.829 $\pm$ 0.030 |  0.668 $\pm$ 0.048 |  0.672 $\pm$ 0.010 |  0.578 $\pm$ 0.011
Absence-indication |  0.882 $\pm$ 0.004 |  0.742 $\pm$ 0.003 |  0.704 $\pm$ 0.010 |  0.589 $\pm$ 0.004
Table 4: Comparison of metrics when calculating attentions between different
views in LANET.

#### The contribution of different feature embeddings to model performance.

The input to our LANET model consists of a label, time, amount, and ID that we
embed to input to the self-attention layer. All these parts include
information on the events preceding the considered one. We want to understand
the contribution of each feature embedding to the model prediction. So, in
each experiment, we drop either amount, time, and ID one at a time. A
significant quality drop after removal indicates the importance of the model
part for the prediction. The results of such an experiment are in Table 5.

Inputs |  Micro-AUC$\uparrow$ |  Macro-AUC$\uparrow$ |  Micro-F1$\uparrow$ |  Macro-F1$\uparrow$
---|---|---|---|---
All |  0.881 $\pm$ 0.007 |  0.737 $\pm$ 0.017 |  0.704 $\pm$ 0.018 |  0.591 $\pm$ 0.003
No amount |  0.825 $\pm$ 0.051 |  0.705 $\pm$ 0.039 |  0.698 $\pm$ 0.010 |  0.574 $\pm$ 0.014
No time |  0.869 $\pm$ 0.007 |  0.721 $\pm$ 0.020 |  0.698 $\pm$ 0.004 |  0.590 $\pm$ 0.003
No ID |  0.880 $\pm$ 0.003 |  0.732 $\pm$ 0.006 |  0.703 $\pm$ 0.022 |  0.588 $\pm$ 0.004
Table 5: How different types of embeddings in the architecture affect LANET.

A substantial drop in quality happens when we drop the amount embeddings. We
attribute this effect to the fact that amounts are essential for quantifying
the scale of interdependencies. Time embeddings are helpful and make models
capture label connections better. IDs contribute little to the performance of
the model. It seems that it is enough to look at other types of available
information.

#### Dependence of model quality on $\tau$.

The model considers only $\tau$ last timestamps for prediction. Thus, a
natural question arises about how the model quality depends on the number of
timestamps used for prediction. The dependence of micro-AUC metric on $\tau$
is presented in 3(a).

(a) The dependence of micro-AUC on the $\tau$ parameter.

(b) The dependence of micro-AUC on the dimension of embeddings.

Figure 3: The dependence of metrics on model parameters for Demand dataset.

The experiment supports the evidence that the more information we provide to
our model, the better prediction we obtain for the considered timestamp.
However, we should find the balance between the desired quality and the
computational cost since the number of involved labels increases with boosting
parameter value. Therefore, for the optimality of the solution and the speed
of the method, we selected the parameter $\tau=3$.

#### Dependence of the model performance on embedding size.

As we are working with learnable embeddings for label, amount, time, and ID
features, we need to set the dimensionality of their representations. For each
feature, we construct embeddings of the same size. The dimensionality of the
embedding vector defines the capacity of information that a model should
capture and the overall number of learnable parameters. The results are in
3(b). We can observe that the model can not effectively learn large-dimension
representations.

#### Attention matrix for vector representations of labels

Additionally, we constructed an attention matrix for the data labels. This
matrix reflects label importances for making predictions. The model pays
attention to various sequence objects with different weights. Figure 4 shows
three attention matrices for three datasets: Order, Transactions, and Liquor.
The size of the matrix is the number of unique labels in the dataset. We
notice that the attention matrix has a clear dominance of the certain labels
encountered in the sequence over those that are not in it, which is clearly
expressed through the weights. Looking deeper, we see that small-scale
variations in attention describe the connection between particular event
types.

Figure 4: Attention matrix studies for the Orders, Transactions and Liquor
datasets from left to right. They show the representation importances.

##  5 Conclusions and discussion

We consider the sequential multi-label classification problem: given previous
multi-label vectors, what a multi-label vector will be next. For this problem,
we propose the LANET architecture. It explores a connection between past
multi-label vectors in a specific architecture. In particular, LANET uses the
self-attention between labels to encourage capturing historical label
interdependencies. In our pipeline, we can apply any model architecture based
on any similar self-attention. So, it can easily absorb attention-related
innovations making it possible to add interpretability [24] or efficiency [30]
to a model.

LANET demonstrates the best metrics on the five considered datasets. In
particular, it improves the micro-AUC metric from $0.7501$ to $0.9536$ on one
of the datasets. However, the current approach is limited to the case where no
long-term memory is required for prediction, as we aggregate only several
recent timestamps and general embeddings of a particular sequence. Further
studies are welcomed to fill this gap.

It is also interesting to note that the proposed approach naturally fits into
the paradigm of self-supervised learning and directly outputs an embedding of
a client in general and at a particular time. The intricacy of the self-
supervised problem is easy to tune by changing the used time horizon or
including masking during training [12, 14]. By aggregating labels collected
during a day, we also make learning event sequences models more stable,
avoiding the over-complicated problem of predicting the next event type [3]
and challenging likelihood optimization for neural temporal points processes
[25].

## References

  * [1] https://pytorch.org/docs/stable/generated/torch.nn.TransformerEncoderLayer.html
  * [2] Prediction of client gender on card transactions. https://www.kaggle.com/c/python-and-analyze-data-final-project/data (2019)
  * [3] Babaev, D., Ovsov, N., Kireev, I., Ivanova, M., Gusev, G., Nazarov, I., Tuzhilin, A.: Coles: Contrastive learning for event sequences with self-supervision. In: Proceedings of the 2022 International Conference on Management of Data. pp. 1190–1199 (2022)
  * [4] Company, C.: Predict Future Sales. https://www.kaggle.com/c/competitive-data-science-predict-future-sales/data (2018)
  * [5] Du, J., Chen, Q., Peng, Y., Xiang, Y., Tao, C., Lu, Z.: Ml-net: multi-label classification of biomedical texts with deep neural networks. Journal of the American Medical Informatics Association 26(11), 1279–1285 (2019)
  * [6] Durand, T., Mehrasa, N., Mori, G.: Learning a deep convnet for multi-label classification with partial labels. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 647–657 (2019)
  * [7] FELIXZHAO: Forecasts for Product Demand. https://www.kaggle.com/datasets/felixzhao/productdemandforecasting (2017)
  * [8] Fursov, I., Morozov, M., Kaploukhaya, N., Kovtun, E., Rivera-Castro, R., Gusev, G., Babaev, D., Kireev, I., Zaytsev, A., Burnaev, E.: Adversarial attacks on deep models for financial transaction records. In: Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining. pp. 2868–2878 (2021)
  * [9] Hang, J.Y., Zhang, M.L.: Collaborative learning of label semantics and deep label-specific features for multi-label classification. IEEE Transactions on Pattern Analysis and Machine Intelligence (2021)
  * [10] Hartvigsen, T., Sen, C., Kong, X., Rundensteiner, E.: Recurrent halting chain for early multi-label classification. In: Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. pp. 1382–1392 (2020)
  * [11] Hawkes, A.G.: Spectra of some self-exciting and mutually exciting point processes. Biometrika 58(1), 83–90 (1971)
  * [12] He, K., Chen, X., Xie, S., Li, Y., Dollár, P., Girshick, R.: Masked autoencoders are scalable vision learners. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16000–16009 (2022)
  * [13] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural computation 9(8), 1735–1780 (1997)
  * [14] Kenton, J.D.M.W.C., Toutanova, L.K.: Bert: Pre-training of deep bidirectional transformers for language understanding. In: Proceedings of NAACL-HLT. pp. 4171–4186 (2019)
  * [15] Lanchantin, J., Wang, T., Ordonez, V., Qi, Y.: General multi-label image classification with transformers. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). pp. 16478–16488 (June 2021)
  * [16] Li, J., Wang, Y., McAuley, J.: Time interval aware self-attention for sequential recommendation. In: Proceedings of the 13th international conference on web search and data mining. pp. 322–330 (2020)
  * [17] Li, Y., Song, Y., Luo, J.: Improving pairwise ranking for multi-label image classification. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (July 2017)
  * [18] Liu, M., Ren, S., Ma, S., Jiao, J., Chen, Y., Wang, Z., Song, W.: Gated transformer networks for multivariate time series classification. arXiv preprint arXiv:2103.14438 (2021)
  * [19] Liu, W., Wang, H., Shen, X., Tsang, I.: The emerging trends of multi-label learning. IEEE transactions on pattern analysis and machine intelligence (2021)
  * [20] Mei, H., Yang, C., Eisner, J.: Transformer embeddings of irregularly spaced events and their participants. In: International Conference on Learning Representations (2021)
  * [21] Menon, A.K., Rawat, A.S., Reddi, S., Kumar, S.: Multilabel reductions: what is my loss optimising? Advances in Neural Information Processing Systems 32 (2019)
  * [22] Pal, A., Selvakumar, M., Sankarasubbu, M.: Multi-label text classification using attention-based graph neural network. arXiv preprint arXiv:2003.11644 (2020)
  * [23] Quadrana, M., Cremonesi, P., Jannach, D.: Sequence-aware recommender systems. ACM Computing Surveys (CSUR) 51(4), 1–36 (2018)
  * [24] Serrano, S., Smith, N.A.: Is attention interpretable? In: Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics. pp. 2931–2951 (2019)
  * [25] Shchur, O., Türkmen, A.C., Januschowski, T., Günnemann, S.: Neural temporal point processes: A review. arXiv preprint arXiv:2104.03528 (2021)
  * [26] Shou, X., Gao, T., Subramaniam, S., Bhattacharjya, D., Bennett, K.: Concurrent multi-label prediction in event streams. In: AAAI Conference on Artificial Intelligence (2023)
  * [27] T, B.: Time Series Analysis of Iowa Liquor Sales. https://www.kaggle.com/code/roudradas/time-series-analysis-of-iowa-liquor-sales/data (2020)
  * [28] Tang, J., Wang, K.: Personalized top-n sequential recommendation via convolutional sequence embedding. In: Proceedings of the eleventh ACM international conference on web search and data mining. pp. 565–573 (2018)
  * [29] Tarekegn, A.N., Giacobini, M., Michalak, K.: A review of methods for imbalanced multi-label classification. Pattern Recognition 118, 107965 (2021)
  * [30] Tay, Y., Dehghani, M., Bahri, D., Metzler, D.: Efficient transformers: A survey. ACM Computing Surveys 55(6), 1–28 (2022)
  * [31] Tsai, C.P., Lee, H.Y.: Order-free learning alleviating exposure bias in multi-label classification. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 34, pp. 6038–6045 (2020)
  * [32] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [33] Wang, S., Hu, L., Cao, L.: Perceiving the next choice with comprehensive transaction embeddings for online recommendation. In: Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2017, Skopje, Macedonia, September 18–22, 2017, Proceedings, Part II 17. pp. 285–302. Springer (2017)
  * [34] Wang, S., Hu, L., Wang, Y., Cao, L., Sheng, Q.Z., Orgun, M.: Sequential recommender systems: challenges, progress and prospects. In: 28th International Joint Conference on Artificial Intelligence, IJCAI 2019. pp. 6332–6338. International Joint Conferences on Artificial Intelligence (2019)
  * [35] Wu, C.Y., Ahmed, A., Beutel, A., Smola, A.J., Jing, H.: Recurrent recommender networks. In: Proceedings of the tenth ACM international conference on web search and data mining. pp. 495–503 (2017)
  * [36] Wu, X.Z., Zhou, Z.H.: A unified view of multi-label performance measures. In: International Conference on Machine Learning. pp. 3780–3788. PMLR (2017)
  * [37] Xiao, L., Huang, X., Chen, B., Jing, L.: Label-specific document representation for multi-label text classification. In: Proceedings of the 2019 conference on empirical methods in natural language processing and the 9th international joint conference on natural language processing (EMNLP-IJCNLP). pp. 466–475 (2019)
  * [38] Yazici, V.O., Gonzalez-Garcia, A., Ramisa, A., Twardowski, B., Weijer, J.v.d.: Orderless recurrent models for multi-label classification. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 13440–13449 (2020)
  * [39] Yeh, C.K., Wu, W.C., Ko, W.J., Wang, Y.C.F.: Learning deep latent space for multi-label classification. In: Thirty-first AAAI conference on artificial intelligence (2017)
  * [40] Ying, H., Zhuang, F., Zhang, F., Liu, Y., Xu, G., Xie, X., Xiong, H., Wu, J.: Sequential recommender system based on hierarchical attention network. In: IJCAI International Joint Conference on Artificial Intelligence (2018)
  * [41] Zhang, M.L., Zhou, Z.H.: A review on multi-label learning algorithms. IEEE transactions on knowledge and data engineering 26(8), 1819–1837 (2013)
  * [42] Zhang, S., Tay, Y., Yao, L., Sun, A., An, J.: Next item recommendation with self-attentive metric learning. In: Thirty-Third AAAI Conference on Artificial Intelligence. vol. 9 (2019)
  * [43] Zhang, W., Jha, D.K., Laftchiev, E., Nikovski, D.: Multi-label prediction in time series data using deep neural networks. arXiv preprint arXiv:2001.10098 (2020)
  * [44] Zhuzhel, V., Grabar, V., Boeva, G., Zabolotnyi, A., Stepikin, A., Zholobov, V., Ivanova, M., Orlov, M., Kireev, I., Burnaev, E., Rivera-Castro, R., Zaytsev, A.: Continuous-time convolutions model of event sequences (2023)