# A variational Monte Carlo approach to the SU(4) spin-orbital model on the
triangular lattice

Chun Zhang  Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  Hui-Ke Jin  Technical University of Munich, TUM School of Natural
Sciences, Physics Department, 85748 Garching, Germany  Yi Zhou
<EMAIL_ADDRESS>Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China  Kavli Institute for Theoretical Sciences, University of Chinese
Academy of Sciences, Beijing 100190, China  CAS Center for Excellence in
Topological Quantum Computation, University of Chinese Academy of Sciences,
Beijing 100190, China

###### Abstract

Previous investigations have suggested that the simplest spin-orbital model on
the simplest frustrated lattice can host a nematic quantum spin-orbital liquid
state. Namely, the orbital degeneracy of the SU(4) Kugel-Khomskii (KK) model
tends to enhance quantum fluctuations and stabilize a quantum spin-orbital
liquid exhibiting stripy features on the triangular lattice, as revealed by
the state-of-the-art method of the density matrix renormalization group
boosted by Gutzwiller projected wave functions. In this work, using the
variational quantum Monte Carlo method, we have studied several spin-orbital
liquid states, including a uniform $\pi$ flux state, three stripy states, and
a plaquette state, on the $L\times{}L$ torus up to $L=24$. It turns out that
one of these stripy states, called the “stripe-II” state, is energetically
favored. This ground state breaks the $C_{6}$ symmetry of the lattice,
resulting in a reduced $C_{2}$ symmetry and doubled unit cells, while
preserving the SU(4) spin-orbital rotation symmetry. Such a nematic quantum
spin-orbital liquid state can be characterized by a parton Fermi surface (FS)
consisting of open orbits in the Brillouin zone, in contrast to the circular
FS of the uniform $\pi$-flux state.

##  I Introduction

The search for quantum spin liquids (QSLs) is one of the central issues in
modern condensed matter physics [1, 2, 3, 4, 5, 6, 7]. In the 1970s, P. W.
Anderson proposed a disordered ground state for the spin-1/2 antiferromagnetic
(AFM) Heisenberg model on the triangular lattice [1], namely the simplest spin
model on the simplest geometrically frustrated lattice, which revealed the age
of QSLs and resonating valence-bond (RVB) states. However, after extensive
studies, the academic community has reached a consensus that the spin-1/2 AFM
Heisenberg model on the triangular lattice hosts a 120∘ magnetically ordered
ground state, contrary to Anderson’s original proposal [8, 9, 10]. It means
that the smallest spin quantum $S=1/2$, which gives rise to the strongest
quantum spin fluctuations, is still insufficient to completely suppress the
classical magnetic order on the triangular lattice, if only the nearest
neighbor (NN) Heisenberg spin interaction $\vec{S}_{i}\cdot\vec{S}_{j}$ is
considered. Therefore, one has to turn to other ways of enhancing quantum
fluctuations to realize a QSL ground state on the simplest geometrically
frustrated lattice.

In addition to (i) the small spin quanta $S$ and (ii) geometric frustration,
there are several other ways to enhance quantum spin fluctuations [11]: (iii)
competing spin interactions involving NN anisotropic spin couplings (e.g., the
anisotropic Ising couplings in the Kitaev honeycomb model [12]) and isotropic
(or anisotropic) spin couplings on longer bonds; (iv) charge fluctuations near
a Mott transition leading to multiple-spin interactions; and (v) additional
degeneracy due to orbital degrees of freedom. Thus, QSL states can be achieved
by stronger geometric frustration [13, 14, 15, 16] or by increasing
competition among mutual interactions [17, 18, 19, 20, 21, 22, 23, 24]. In
particular, orbital degeneracy [i.e. the path (v) mentioned above] can
significantly enhance quantum fluctuations [25]. For instance, in some
transition metal oxides, spin and orbital degrees of freedom may play a
symmetric role, resulting in an enlarged SU(4) rather than the
SU(2)$\times$SU(2) spin-orbital rotational symmetry. The SU(4) Kugel-Khomskii
(KK) model is the minimal model used to describe such spin-orbital materials
[26, 27, 28, 29, 30, 31, 32, 33], in which the SU(4) symmetry can amplify
quantum fluctuations and potentially stabilize a spin-orbital liquid ground
state.

Much theoretical effort has been devoted to the study of quantum phases in
SU(4) quantum magnets. As an example, the one-dimensional SU(4) KK model is
found to be integrable and to possess gapless excitations [34, 35]. These
gapless excitations are well characterized by the low-energy effective theory,
i.e. the SU(4)1 Wess-Zumino-Witten (WZW) conformal field theory (CFT) [36, 37,
38]. The ground state of the SU(4) KK model on the two-leg ladder is also well
known, which breaks the translational symmetry and forms an SU(4) singlet
plaquette [39, 40]. This SU(4) singlet plaquette state can even be an exact
ground state on the ladder, as long as the extra interactions have been
properly chosen and added to the model Hamiltonian [41]. Unlike one-
dimensional (1D) and ladder models, two-dimensional (2D) models are less
clear, despite extensive research. For instance, several complementary methods
indicated that the ground state of the SU(4) KK model on the honeycomb lattice
is a Dirac-type spin-orbital liquid [42, 43, 44]. For the square lattice,
several candidate ground states have been proposed for the SU(4) KK model,
such as the plaquette ordered state [27, 45], the $\mathbf{Z}_{2}$ spin-
orbital liquid [46], and the dimerized and SU(4) symmetry-breaking state [47].

Especially, the ground state of the SU(4) KK model on the triangular lattice
has been debated until very recently. In the early days, exact diagonalization
(up to 16 sites) and variational approach (up to 64 sites) suggested an
SU(4)-singlet plaquette liquid ground state on which simple types of long-
range correlations are suppressed [48]. Recently, the combination study of
density matrix renormalization group (DMRG) calculation on cylinders (with a
circumference up to $L_{y}=4$) and field theory analysis proposed a gapless
liquid state with an emergent parton Fermi surface (FS) in the 2D limit [49].
On the contrary, another DMRG study [50] and a self-consistent mean-field
theory [51] prefer a stripe-ordered ground state. Soon thereafter, using a
state-of-the-art method, DMRG boosted by Gutzwiller projected wave functions
[52], two of the authors and their collaborators revisited this model [53].
The latest work revealed a nematic quantum spin-orbital liquid state, a
critical stripy state that preserves SU(4) symmetry but breaks translational
symmetry by doubling the unit cell along one of two primitive vectors. It was
shown that the central charge of each stripe is $c=3$, which is in
quantitative agreement with the SU(4)1 WZW CFT. Moreover, it was found that
the DMRG-obtained state can be well described by a “single” Gutzwiller
projected wave function with an emergent parton FS that consists of open
orbits in the reciprocal space and indicates the nematicity. Such a “unified”
picture captures all the essential physics in both quasi-1D cylinders (with
circumferences up to $L_{y}=8$) and the 2D limit, although the former exhibits
a strong finite-size effect and even-odd discrepancy [53]. Note that the newly
proposed parton FS consists of open orbits in reciprocal space, which must
undergo a Lifshitz transition from the originally proposed $\pi$-flux state,
which has a closed parton FS [49].

In this work, we re-examine the SU(4) KK model on the triangular lattice by
using the variational Monte Carlo (VMC) approach. The advantage of the VMC
method is to reduce the impact of finite-size effects in 2D systems and to
impose periodic boundary conditions (PBCs) in both directions. Thus, the
previously proposed Gutzwiller projected wave function [53] can be carefully
verified at large tori, up to a lattice size of $24\times{}24$. To this end,
we have proposed five types of spin-orbital liquid states: a uniform
$\pi$-flux state, three stripy states, and a plaquette state. After extensive
numerical efforts, we have gathered strong evidence pointing to a critical
stripy ground state, as revealed by the Gutzwiller-boosted DMRG [53].

The rest of this paper is organized as follows. We revisit the SU(4) KK model
in Section II to make it self-contained. Then five fermionic parton mean-field
ansatzes have been proposed in Section III. In Section IV, we search for the
ground state with the help of VMC and the stochastic reconfiguration method.
Section V is devoted to summary.

##  II Model and symmetries

The SU(4) KK model Hamiltonian on the triangular lattice is defined as

$H=\frac{1}{2}\sum_{\langle
i,j\rangle}\left(4\vec{S}_{i}\cdot\vec{S}_{j}+1\right)\left(4\vec{T}_{i}\cdot\vec{T}_{j}+1\right),$
(1a) where $\langle i,j\rangle$ denotes an NN bond and $\vec{S}_{i}$
($\vec{T}_{i}$) is the $S=1/2$ spin (orbital) vector at site $i$. These spin
(orbital) vectors can be represented by introducing
$\vec{S}_{i}=\frac{1}{2}\\{\sigma^{x}_{i},\sigma^{y}_{i},\sigma^{z}_{i}\\}$
($\vec{T}_{i}=\frac{1}{2}\\{\tau^{x}_{i},\tau^{y}_{i},\tau^{z}_{i}\\}$), where
$\sigma^{x,y,z}$ ($\tau^{x,y,z}$) are the standard Pauli matrices acting on
the two-fold spin (orbital) indices. It is easy to verify that the Hamiltonian
$H$ in Eq. (1) commutes with all the 15 generators of the SU(4) Lie group,
$\lambda^{1},\lambda^{2},\cdots,\lambda^{15}$, that are linear combinations of
$\\{\vec{\sigma},\vec{\tau},\vec{\sigma}\otimes{}\vec{\tau}\\}$. Therefore the
SU(4) symmetry is respected by the model Hamiltonian.

The SU(4) symmetry becomes more transparent from the point of view that the
$2\times{}2=4$ spin-orbital degrees of freedom at each site can be treated as
a pseudospin. To demonstrate this, we introduce a 15-dimensional vector
$\vec{\lambda}_{i}\equiv\\{\lambda_{i}^{1},\lambda_{i}^{2},\cdots,\lambda_{i}^{15}\\}$
whose components are 15 SU(4) generators. Then the model Hamiltonian can be
rewritten in terms of $\vec{\lambda}_{i}$ as

$H=\sum_{\langle
i,j\rangle}\left(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j}+\frac{1}{2}\right).$
(1b)

Note that these generators, $\lambda_{i}^{\mu}$’s, have been normalized with
$\mbox{tr}(\lambda_{i}^{\mu}\lambda_{j}^{\nu})=2\delta_{ij}\delta_{\mu\nu}$
[54]. So the Hamiltonian in Eq. (1) can be interpreted as an SU(4) AFM
Heisenberg model, which is invariant under the global SU(4) pseudospin-
rotation symmetry.

In addition to the SU(4) symmetry, this model is also symmetric with respect
to the spatial symmetries of the triangular lattice, including the lattice
translation symmetries $T_{1,2}$, the mirror symmetry $M$, and the sixfold
rotation symmetry $C_{6}$. Furthermore, as a spin-orbital system, this model
naturally preserves the time-reversal symmetry.

##  III Parton construction and mean-field ansatzes

In this section, we shall first use the Gutzwiller projection to construct a
trial wave function for solving the Hamiltonian $H$ in Eq. (1). Given a
specific variational wave function $|\Psi_{trial}\rangle$, one can calculate
the ground-state energy $E$,

$E=\frac{\langle\Psi_{trial}{}|H|\Psi_{trial}\rangle}{\langle\Psi_{trial}|\Psi_{trial}\rangle},$

by using the standard Monte Carlo method. This trial wave function can then be
optimized by minimizing the energy $E$.

To construct such a trial Gutzwiller projected state, we introduce four
flavors of fermionic partons at each site $i$, $f^{\dagger}_{i\alpha}$ for
$\alpha=1,2,3,4$ ($f_{i\alpha}$ the annihilation operators), where $\alpha$
denotes the four pseudospin (spin-orbital) degrees of freedom mentioned above.
The SU(4) pseudospin operators can now be expressed in terms of these
fermionic partons as

$\lambda^{\mu}_{i}=f_{i}^{\dagger}\lambda^{\mu}f_{i},\quad{\rm~{}for~{}}\mu=1,\cdots,15,$
(2)

where $f_{i}\equiv(f_{i1},f_{i2},f_{i3},f_{i4})^{T}$ is a four-component
vector. Instead of $4$ spin-orbital states, these four flavors of fermions
will generate $2^{4}=16$ states in Fock space at each local site. To restore
the physical Hilbert space, a single occupancy constraint of
$\sum_{\alpha=1}^{4}f_{i\alpha}^{\dagger}f_{\alpha}=1$ was imposed on each
site, e.g. the Gutzwiller projection was implemented.

In this parton representation in Eq. (2), the inner product of two pseudospins
$\vec{\lambda}_{i}\cdot\vec{\lambda}_{j}$ can be written in terms of four-
fermion interactions:

$\begin{split}\vec{\lambda}_{i}\cdot\vec{\lambda}_{j}=&-2\sum_{\alpha,\beta=1}^{4}\left({}f_{i\alpha}^{{\dagger}}f_{j\alpha}f_{j\beta}^{{\dagger}}f_{i\beta}+\frac{1}{4}{}f_{i\alpha}^{{\dagger}}f_{i\alpha}f_{j\beta}^{{\dagger}}f_{j\beta}\right)\\\
&+\sum_{\alpha=1}^{4}\left(f^{\dagger}_{i\alpha}f_{i\alpha}+f^{\dagger}_{j\alpha}f_{j\alpha}\right)\end{split}.$
(3)

Here the identity
$\vec{\lambda}_{ab}\cdot\vec{\lambda}_{cd}=2\delta_{ad}\delta_{bc}-\frac{1}{2}\delta_{ab}\delta_{cd}$
has been used. Thus, the model Hamiltonian $H$ in Eq. (1) can be recast as
follows:

$H=\sum_{\langle{}ij\rangle}\left(-\sum_{\alpha,\beta=1}^{4}{}f_{i\alpha}^{{\dagger}}f_{j\alpha}f_{j\beta}^{{\dagger}}f_{i\beta}+1\right).$
(4)

Note that the single-occupancy condition was imposed in the derivation of the
last equation. The four-fermion term can be decoupled by introducing the mean-
field parameters on NN bonds $\langle{}ij\rangle$ as

$\chi_{ij}=\chi_{ji}^{\ast}=\sum_{\alpha=1}^{4}\langle
f_{i\alpha}^{\dagger}f_{j\alpha}\rangle.$ (5)

So we get a quadratic mean-field Hamiltonian

$H_{MF}=\sum_{\langle{}ij\rangle}\sum_{\alpha=1}^{4}\chi_{ij}f_{j\alpha}^{\dagger}f_{i\alpha}+h.c.,$
(6)

which preserves the SU(4) symmetry and can be easily diagonalized now.

In fact, $H_{MF}$ in Eq. (6) can be treated as a “variational Hamiltonian” or
an “effective Hamiltonian” on which $\\{\chi_{ij}\\}$ is a set of variational
parameters to be determined. Taking into account the single occupancy
constraint at the mean, i.e.
$\sum_{\alpha=1}^{4}\langle{}f_{i\alpha}^{\dagger}f_{i\alpha}\rangle=1$, one
can find a ground state $|\Psi_{MF}\rangle$ for the mean-field Hamiltonian
$H_{MF}$ at $1/4$ filling. Then the trial wave function $|\Psi_{trial}\rangle$
can be constructed by performing the Gutzwiller projection on
$|\Psi_{MF}\rangle$:

$|\Psi_{trial}\rangle=P_{G}|\Psi_{MF}\rangle.$ (7)

Here $P_{G}$ is the Gutzwiller projector, which removes all non-single-
occupied components to locally enforce the single-occupancy constraint.

The set of parameters, $\\{\chi_{ij}\\}$, together with the form of $H_{MF}$
in Eq. (6), is called “mean-field Ansatz”. As long as a mean-field Ansatz is
given, one can evaluate the energy $E(\\{\chi_{ij}\\})$ by using the standard
VMC method, and minimize $E(\\{\chi_{ij}\\})$ to obtain an optimal set of
parameters $\\{\chi_{ij}\\}$. Then, all the correlation functions can be
computed with such an optimized wave function. It is worth mentioning that the
mean-field parameters $\chi_{ij}$ can be determined up to an overall positive
factor via the ground state energy optimization. In the remaining part of this
section, we shall discuss several typical examples of SU(4)-symmetric mean-
field Ansatzes which will be utilized to construct trial wave functions.

Figure 1: For a $2\times{}2$ VBC on a triangular lattice, the unit cell has
been enlarged to include four sublattices denoted by A, B, C, and D,
respectively. (a) The XC torus geometry and the periodic boundary condition
(PBC) are shown by dashed lines. (b) The 12 types of NN bonds in the enlarged
unit cell.

Uniform $\pi$-flux state: Two of the simplest mean-field ansatzes are uniform
zero-flux or $\pi$-flux states. Since the latter always has a lower energy
than the former, we will skip the higher-energy zero-flux state and discuss
only the uniform $\pi$-flux state. The uniform $\pi$-flux state is governed by
the mean-field Hamiltonian as follows,

$H_{\Delta=\pi}=-t\sum_{\alpha=1}^{4}\sum_{\langle
ij\rangle}(f_{i\alpha}^{\dagger}f_{j\alpha}+h.c.),~{}$ (8)

which is obtained by setting $\chi_{ij}=-t\,\,(t>0)$ on all the NN bonds in
Eq. (6). In each elementary triangle, $\chi_{ij}$ will pick up a $\pi$-flux
that gives a name to the $\pi$-flux state proposed in Ref. [49]. The parton FS
at $1/4$ filling is plotted in Fig. 2(a), which is almost a circle.

As an instability of the uniform $\pi$-flux state, a “stripy state” has
recently been proposed [53]. The stripy state breaks the lattice translational
symmetry by doubling the unit cell along a certain lattice vector direction.
In this work, we consider a more generic class of valence-bond-crystal (VBC)
states. These VBC states have an enlarged $2\times{}2$ unit cell, as
illustrated in Fig. 1(a). One can obtain stripy states or the uniform
$\pi$-flux state by restoring the lattice translational symmetry along one or
two lattice vector directions.

Figure 2: The parton energy contour plot in the first (unfolded) Brillouin
zone (BZ) for (a) the uniform $\pi$-flux state [Eq. (8)], (b) the stripe-I
state [$\delta_{1}=0.15$ in Eq. (12)], (c) the stripe-II state
[$\delta_{1}=0.15$ and $\epsilon=0$ in Eq. (13)], (d) the stripe-III state
[$\delta_{1}=0.15$ in Eq. (14)], and (e) the plaquette state
[$\delta_{1}=\delta_{2}=0.22$ in Eq. (15)], where the parameter $t=1$ be set.
The dashed lines in (b)-(e) enclose the folded BZ. At $1/4$ filling, (a) the
uniform $\pi$-flux state has a closed Fermi surface (FS), (b)-(d) stripe
states allow an open FS in the BZ, while (e) the plaquette state is a parton
band insulator.

Generic $2\times{}2$ VBC states: Since the coordination number is $z=6$ on a
triangular lattice. Each unit cell in a $2\times{}2$ VBC state has
$2\times{}2\times{}6/2=12$ NN bonds, as labeled by $n=1,2,\cdots,12$ in Fig.
1(b). Thus, a generic VBC state can be depicted by 12 mean-field parameters
$\chi_{ij}$ within one unit cell as

$\chi_{ij}=-t_{n},$ (9)

such that the corresponding mean-field Hamiltonian takes the form of

$H_{VBC}=-\sum_{\alpha=1}^{4}\sum_{n=1}^{12}\sum_{\langle
ij\rangle_{n}}t_{n}\left(f_{i\alpha}^{\dagger}f_{j\alpha}+h.c.\right),~{}$
(10)

where $\langle{}ij\rangle_{n}$ denotes an type-$n$ NN bond as shown in Fig.
1(b). As mentioned above, stripy states can be obtained by restoring the
lattice translational symmetry along either the $x$\- or the $y$-direction.

Table 1: Five types of mean-field ansatzes. Here $t_{n}>0$ on all the NN bonds $\langle{ij}\rangle_{n=1,\cdots,12}$, suggesting that $\pi$-flux states are energetically favorable. $T_{x}(T_{y})$ shifts the system along the $x(y)$\- direction by a lattice constant.  | schematic bond strength | mean-field ansatz | unit cell | symmetry
---|---|---|---|---
(a) uniform $\pi$-flux state |  |  $\chi_{ij}=-t$ see Eq. (8) | $1\times{}1$ |  $T_{x},T_{y}$ $D_{6}$
(b) stripe-I state |  |  $t_{n+6}=t_{n}|_{(n=1,\cdots,6)}$ $t_{2}=t_{4}=t_{5}$ see Eq. (12) | $2\times{}1$ |  $T_{x}^{2},T_{y}$ $C_{2}$
(c) stripe-II state |  |  $t_{n+6}=t_{n}|_{(n=1,\cdots,6)}$ $t_{2}=t_{5}$ see Eq. (13) | $2\times{}1$ |  $T_{x}^{2},T_{y}$ $C_{2}$
(d) stripe-III state |  |  $t_{n+6}=t_{n}|_{(n=1,\cdots,6)}$ $t_{1}=t_{2}$ $t_{4}=t_{5}$ see Eq. (14) | $2\times{}1$ |  $T_{x}^{2},T_{y}$ $C_{2}$
(e) plaquette state |  |  $t_{1}=t_{2}=t_{3}$ $t_{4}=t_{5}=t_{6}$ $t_{7}=t_{8}=t_{9}$ $t_{10}=t_{11}=t_{12}$ $t_{1}+t_{10}=t_{4}+t_{7}$ see Eq. (15) | $2\times{}2$ | $T_{x}^{2},T_{y}^{2}$

$C_{2}$\- stripe states: From a symmetry point of view, a stripy state will
break the lattice translational symmetry in one direction while preserving the
other. The $D_{6}$ rotational symmetry of the lattice will be broken
simultaneously. The resulting stripy state may or may not remain a reduced
$C_{2}$ rotational symmetry of the original $D_{6}$ symmetry. In particular,
we are interested in the stripy states exhibiting the $C_{2}$ rotational
symmetry, dubbed the “$C_{2}$\- stripe states”.

Without loss of generality, we focus on stripy states that only break the
translational symmetry along the $x$-direction (associated with symmetry
generator $T_{x}$). The unit cell is of $2\times{}1$ and consists of 6 NN
bonds. Therefore, there are only two sublattices, e.g., $\mbox{A}=\mbox{C}$
and $\mbox{B}=\mbox{D}$ in Fig. 1(b). The 12 values of $\chi_{ij}$, say,
$t_{n}$ ($n=1,2,\cdots,12$) defined according to the labeling in Fig. 1(b) are
reduced to 6 values in such a $2\times{}1$ VBC state. Namely, the constraint
$t_{n}=t_{n+6}$ ($n=1,2,3,4,5,6$) should be imposed on Eqs. (9) and (10).
Introduce a lattice site labeling scheme,

$i=(x_{i},y_{i})=x_{i}\hat{x}+y_{i}\hat{y},$

we find that for a $2\times{}1$ stripy state, $\chi_{ij}$ must satisfy on NN
bonds:

$\chi_{i,i+\hat{x}}=-t_{1(4)},\,\,\,\chi_{i,i+\hat{y}}=-t_{2(5)},\,\,\,\chi_{i,i+\hat{x}-\hat{y}}=-t_{3(6)},\,\,\,i\in\mbox{A
(B)}.$

Now we go further to consider $C_{2}$\- stripe states that respect the reduced
$C_{2}$ rotational symmetry around the $\hat{x}-\frac{1}{2}\hat{y}$ axis,
giving rise to $t_{1}=t_{3}$ and $t_{4}=t_{6}$. Therefore, the mean-field
ansatz will take a simplified form on the 6 NN bonds in the $2\times{}1$ unit
cell, as

$\chi_{i,i+\hat{x}}=\chi_{i,i+\hat{x}-\hat{y}}=-t_{1(4)},\,\,\,\chi_{i,i+\hat{y}}=-t_{2(5)},\,\,\,i\in\mbox{A
(B)}.$ (11)

Note that due to symmetry there are only 4 different values of $\chi_{ij}$,
i.e. $t_{1}$, $t_{2}$, $t_{4}$ and $t_{5}$. Moreover, such a $C_{2}$-stripe
state is characterized by at most three independent parameters up to a
positive total factor. In the following, we will discuss three types of
$C_{2}$\- stripe states, called “stripe-I, -II, and -III states”, which have
been found by the VMC optimization to be local minima in the full parameter
space for $2\times{}2$ VBC states, $\\{t_{1},t_{2},\cdots,t_{12}\\}$.

Stripe-I state: This type of states can be achieved by setting

$t_{2}=t_{4}=t_{5}$ (12a) in Eq. (11). The mean-field ansatz for such a stripy
state can be simplified as follows: $\chi_{ij}=-t(1+r_{ij}\delta_{1}),$ (12b)
where $r_{ij}=\pm{}1$ is an integer defined by
$r_{ij}=\delta_{x_{i},x_{j}}+(-1)^{x_{i}}(1-\delta_{x_{i},x_{j}}),$ (12c)

$t=(t_{1}+t_{4})/2$ serves as an overall factor,
$\delta_{1}=(t_{4}-t_{1})/(t_{4}+t_{1})$ is the contrast of the stripe, and
$x_{i}\leq{}x_{j}$ is assumed for NN bond $\langle{}ij\rangle$. Note that the
stripe-I state is exactly equivalent to the trial wave function used to
initialize the DMRG calculations in Ref. [53].

Stripe-II state: Imposing the constraint

$t_{2}=t_{5}$ (13a) in Eq. (11) gives rise to the stripe-II state. The mean-
field ansatz for such a state reads
$\chi_{ij}=-t\left[1+(-1)^{x_{i}}\delta_{1}\right](1-\delta_{x_{i},x_{j}})-t(1+\epsilon)\delta_{x_{i},x_{j}},$
(13b)

where $t=(t_{1}+t_{4})/2$ is an overall factor,
$\delta_{1}=(t_{4}-t_{1})/(t_{4}+t_{1})$ reflects the stripe contrast,
$\epsilon=t_{2}/t-1$, and $x_{i}\leq{}x_{j}$ has been assumed for NN bond
$\langle{}ij\rangle$. When $\epsilon=\delta_{1}$, the stripe-II state becomes
a stripe-I state.

Stipe-III state: Such a state can be achieved by setting

$t_{1}=t_{2}\,\mbox{ and }\,t_{4}=t_{5}$ (14a) in Eq. (11). The corresponding
mean-field ansatz reads $\chi_{ij}=-t\left[1+(-1)^{x_{i}}\delta_{1}\right],$
(14b)

where $t=(t_{1}+t_{4})/2$ and $\delta_{1}=(t_{4}-t_{1})/(t_{4}+t_{1})$, and
$x_{i}\leq{}x_{j}$ has been assumed for NN bond $\langle{}ij\rangle$.

Plaquette state: In addition to the uniform $\pi$-flux state and the three
$C_{2}$\- stripe states, we also consider the “plaquette state”, which is
given by “perfect triangle conditions”:

$\begin{split}t_{1}&=t_{2}=t_{3},\\\ t_{4}&=t_{5}=t_{6},\\\
t_{7}&=t_{8}=t_{9},\\\ t_{10}&=t_{11}=t_{12},\end{split}$ (15a) and the
constraint $t_{1}+t_{10}=t_{4}+t_{7}.$ (15b) Thus, the mean-field ansatz for
such a plaquette state takes a three-parameter form of
$\chi_{ij}=-t[1+(-1)^{x_{i}}\delta_{1}+(-1)^{\min\\{y_{i},y_{j}\\}}\delta_{2}].$
(15c) where $x_{i}\leq{}x_{j}$ has been assumed for NN bond
$\langle{}ij\rangle$, and the three parameters $t$, $\delta_{1}$, as well as
$\delta_{2}$ are determined by $\begin{split}t_{1}+t_{10}&=t_{4}+t_{7}=2t,\\\
t_{1}&=t(1+\delta_{1}+\delta_{2}),\\\
t_{4}&=t(1+\delta_{1}-\delta_{2}).\end{split}$ (15d)

It is worth noting that the parameter $t$ for these five states that are
defined in Eqs. (8), (12), (13), (14) and (15) respectively can be unified as
the average of the 12 $t_{n}$’s in Eq. (10),

$t=\frac{1}{12}\sum_{n=1}^{12}t_{n}.$ (16)

Before the end of this section, we summarize five types of typical parton
states in Table 1, including (a) the uniform $\pi$-flux state, (b) the
stripe-I state, (c) the stripe-II state, (d) the stripe-III state, and (e) the
plaquette state. With $1/4$ parton filling, all stripy states allow an open
parton FS when the stripe contrast exceeds a critical value,
$\delta_{1}>\delta_{c}$ (with $\delta_{c}\sim{}0.1$), as shown in Fig. 2.
However, the plaquette state is always a band insulator of fermionic partons.

##  IV VMC calculations and results

We have performed VMC calculations on various $L_{x}\times L_{y}$ XC lattices
up to $L_{x}=L_{y}=24$ [see Fig. 1(a)], on which PBCs for SU(4) pseudospins
have been imposed along both $x$\- and $y$-directions. Thus, the parton
boundary conditions can be either PBCs or antiperiodic boundary conditions
(APBCs). Indeed, we have found that the system always has a lower variational
energy when the parton boundary conditions have been chosen to be APBCs rather
than PBCs along both $x$\- and $y$\- directions.

First, we have searched for local energy minima with generic $2\times{}2$ VBC
states given in Eq. (10), parameterized by the set of
$\\{t_{1},t_{2},\cdots,t_{12}\\}$. To handle such a large parameter space,
traditional methods such as steepest descent or Hessian matrix construction
are computationally expensive and inefficient for optimizing the mean-field
parameters. Therefore, a stochastic reconfiguration (SR) method has been
exploited to optimize the parameters [55, 56]. In our calculations, a dynamic
step length for the variational parameters $t_{n}$, $\zeta$, has been
attempted during the VMC simulations, instead of a constant step length. At
the beginning of optimization, a relatively large step length modulation,
$\zeta\sim 5\times 10^{-2}$, is chosen to avoid trapping $E$ in local minima.
Then we gradually decrease the modulation of $\zeta$ and in the last $10^{2}$
Monte Carlo steps, the step length is refined as $\zeta\sim 10^{-4}$. More
details on the SR method can be found in Appendix A.

By performing extensive VMC computations initialized with a large number of
sets of random $\\{t_{1},\cdots,t_{12}\\}$, we can draw two conclusions as
follows:

###### Conclusion 1

All the energetically favorable states, found at local energy minima,
correspond to a $\pi$-flux configuration, i.e., $t_{n}>0$ for all the
$n=1,\cdots,12$ in Eq. (9).

###### Conclusion 2

Within the numerical error bars, each local energy minimum found by the
12-parameter VMC optimizations can be identified as either the uniform
$\pi$-flux state, one of the three $C_{2}$ stripy states, or a plaquette state
(see Table 1).

Based on these two results, we have performed further VMC calculations on
$L_{x}=L_{y}=L$ tori using the five promising parton mean-field ansatzes
listed in Table 1. Explicitly, in addition to the 12-parameter VMC approach,
we have also applied the parameter constraints defined in Eqs. (12), (13),
(14), and (15) to our VMC calculations. Indeed, we find that all these
simulations converge to energy minima that share the same energy with those
obtained by the 12-parameter VMC optimization, supporting Conclusion 2. More
details can be found in Appendix B.

Figure 3: All the local energy minima can be classified into the five types of states listed in Table 1. (a) The optimized ground-state energy for the five mean-field ansatzes listed in Table 1. The lattice geometry is chosen to be an $L\times L$ XC torus with $L=8,12,16,20,24$. (b) The linear fit of $E=E_{\infty}+\beta{}/L^{3}$ with $L=8,12,16,20,24$. The values of $E_{\infty}$ can be found in Table 2.The energy error bar $(\thicksim 10^{-3})$ is approximately equal to the symbol size. Table 2: The optimized ground-state energy and corresponding mean-field parameters for all these five types of local energy minima found in the 12 variational-parameter optimization. $E_{24}$ is the optimized ground state energy on an XC($24\times{}24$) torus. $E_{\infty}$ is obtained from the linear fit $E=E_{\infty}+\beta{}/L^{3}$ with $L=8,12,16,20,24$. The variational parameters $t_{n}(n=1,2,\cdots,12)$, $\delta_{1}$, $\delta_{2}$ and $\epsilon$ are defined in Eqs. (10), (11), (12), (13), and (14). Here we set the average $t=\frac{1}{12}\sum_{n=1}^{12}t_{n}$ as the unit of $t_{n}$ [see Eq. (16)].  | $E_{24}$ | $E_{\infty}$ | $t_{n}$ | $\delta_{1},\delta_{2},\epsilon$
---|---|---|---|---
(a) uniform $\pi$-flux state | $-2.8263\pm 0.0023$ | $-2.8292\pm 0.0025$ |  $t_{n}|_{(n=1,\cdots,12)}\approx t;$ [see Eq. (8)] |  $/$
(b) stripe-I state | $-2.8397\pm 0.0024$ | $-2.8409\pm 0.0023$ |  $t_{1}/t\approx 0.816,\,{}t_{4}/t\approx 1.092;$ [see Eq. (12)] |  |  $\delta_{1}\approx 0.145$
(c) stripe-II state | $-2.8484\pm 0.0026$ | $-2.8500\pm 0.0026$ |  $t_{1}/t\approx 0.849,\,{}t_{4}/t\approx 1.151,$ $t_{2}/t\approx{}1.0;$ [see Eq. (13)] |  $t_{n+6}\approx{}t_{n}|_{(n=1,\cdots,6)}$ $t=(t_{1}+t_{4})/2$ |  $\delta_{1}\approx 0.151$ $\epsilon\approx{}0$
(d) stripe-III state | $-2.8335\pm 0.0022$ | $-2.8349\pm 0.0028$ |  $t_{1}/t\approx 0.854,\,{}t_{4}/t\approx 1.146;$ [see Eq. (14)] |  |  $\delta_{1}\approx 0.146$
(e) plaquette state | $-2.7673\pm 0.0023$ | $-2.7703\pm 0.0026$ |  $t_{1}/t\approx 0.58,\,{}t_{4}\approx t_{7}\approx{}t_{2},\,{}t_{10}/t\approx 1.42,$ $t=(t_{1}+t_{10})/2;$ [see Eq. (15)] |  $\delta_{1}\approx\delta_{2}\approx 0.210$

The main results are summarized in Fig. 3 and Table 2. The optimized ground
state energies have been plotted as a function of the linear size $L$ [see
Fig. 3(a)] as well as its inverse cube $1/L^{3}$ [see Fig. 3(b)] [57, 58, 59].
Note that here the optimizations have been done with the 12 variational
parameters, $t_{1,\cdots,12}$, and all the local energy minima have been
classified into five types of states as listed in Table 1. For all these five
types of local energy minima found in the 12 variational parameter
optimization, the optimized ground state energy and corresponding mean-field
parameters can be found in Table 2.

Using the formula $E(L)=E_{\infty}+\beta/L^{3}$ [57, 58, 59], extrapolation to
the thermodynamic limit $L\to{}\infty$ allows us to draw the third conclusion:

###### Conclusion 3

Among all the five types of ansatzes listed in Table 1, the stripe-II state
has the lowest ground-state energy.

Stability of variational ground states: We will now investigate the energetic
stability of these variational ground states. To do this, we introduce a
Gaussian noise for a given mean-field ansatz $\\{t_{i}\\}$ such that
$t_{n}\rightarrow{}t_{n}+\delta{}t_{n}$, where
${\delta{}t_{n}}\sim{}0.01{}t_{n}$ and the constraint
$\sum_{n}\delta{}t_{n}=0$ has been imposed. We then initialize the VMC
computations with these noisy ansatzes and allow all 12 parameters
$t_{n=1,\cdots,12}$ to be optimized during the VMC computations.

We have observed that the stripe-II state and the plaquette state are stable
against noise in VMC simulations, even when a relatively strong Gaussian noise
$\\{\delta{}t_{i}\\}$ is added. This means that a noisy stripe-II state
converges to a stripe-II state after VMC optimization, and so does a noisy
plaquette state. In contrast, the uniform $\pi$-flux state, the stripe-I
state, and the stripe-III state are inherently unstable to Gaussian noise.
They easily converge to a stripe-II state instead of the original state type
during VMC optimization.

Most stable ground state: We now discuss the most stable ground state, the
stripe-II state defined in eq. (13). First, we illustrate in Fig. 4 the set of
all 12 variational parameters $\\{t_{1},t_{2},\cdots,t_{12}\\}$ for the fully
optimized state on an XC($24\times{}24$) torus. Note that such a state
satisfies the constraints $t_{n+6}=t_{n}|_{n=1,\cdots,6}$ and $t_{2}=t_{5}$
within the error bar, which defines a stripe-II state. Moreover, this stripe-
II state reaches its lowest energy at $t_{1}+t_{4}=2t_{2}$ or $\epsilon=0$ in
eq. (13).

Figure 4: Stipe-II state on an XC($24\times 24$) torus. The lowest energy
state is characterized by 12 variational parameters $t_{n=1,\cdots,12}$, which
are illustrated in a gray color-map. (a) A large unit cell on the XC($24\times
24$) torus for a $2\times{}2$ generic VBC state. The dashed lines indicate
PBCs. (b) A unit cell for the strip-II state. Here $t_{1}+t_{4}=2t_{2}$ or
$\epsilon=0$ [defined in Eq. (13)] is found.

Second, we will focus on the stripe-II state defined in Eq. (13), which is
characterized by two independent parameters $\delta_{1}$ and $\epsilon$.
Motivated by the above result of the stripe-II state during the 12-parameter
optimization, we fix $\epsilon=0$ and explore the ground state energy by
varying $\delta_{1}$ in Eq. (13). The ground state energy
$E(\delta_{1},\epsilon=0)$ was calculated on a $L\times{}L$ torus and plotted
in Fig. 5 as a function of $\delta_{1}$, considering different system sizes of
$L=8,12,16,20,24$. Remarkably, the minimum energy was found to be
$E_{min}=-2.8386\pm{}0.0027$ with $\delta_{1}=0.15$ on a $L=24$ lattice,
showing a strong agreement with the results obtained by the 12-parameter
optimization.

Figure 5: Stipe-II states on an XC($L\times{}L$) torus. Here we have set
$\epsilon=0$ in Eq. (13). The ground state energy $E(\delta_{1},\epsilon=0)$
is plotted as a function of $\delta_{1}$ for $L=8,12,16,20,24$. The energy
error bar $(\thicksim 10^{-3})$ is approximately equal to the symbol size. The
energy minimum is found to be $E_{min}=-2.8386\pm{}0.0027$ with
$\delta_{1}=0.15$ on the $L=24$ lattice.

Finite-size effect: Finally, we would like to examine the finite-size effect
on lattices with different aspect ratios, $L_{x}/L_{y}$. To do this, we fix
$L_{x}$ and vary $L_{y}$, and then perform VMC calculations on the
XC($L_{x}\times{}L_{y}$) tori by the 12-parameter optimization. A number of
typical results are shown in Fig. 6. We find that the stripe-II state remains
the lowest energy state as the aspect ratio $L_{x}/L_{y}$ varies, i.e. our
conclusions remain unaffected by variations in the aspect ratio.

Figure 6: The optimized energy $E(t_{1},\cdots,t_{12})$ for different shapes
of XC($L_{x}\times{}L_{y}$) tori. Here $L_{x}=12$ is fixed and $L_{y}$ is
chosen to be $L_{y}=6,8,12$. Inserts (1), (2), and (3) represent three types
of states in a gray color map. The energy error bar $(\thicksim 10^{-3})$ is
approximately equal to the symbol size.

##  V Summary

In summary, the VMC method has been exploited to study the SU(4) symmetric
spin-orbital model, i.e., the so-called Kugel-Khomskii model, on the
triangular lattice. Begin with a generic $2\times{}2$ valence-bond-crystal
state that is characterized by 12 variational parameters
$\\{t_{1},t_{2},\cdots,t_{12}\\}$, we have performed extensive VMC
calculations and found that (i) all the local energy minima exhibit a
$\pi$-flux configuration, (ii) the optimized local energy minimum must belong
to one of five types of states that include the uniform $\pi$-flux state,
three $C_{2}$\- stripe states (dubbed stripe-I, stripe-II, and stripe-III
states), and the plaquette state, as listed in Table 1, and (ii) the stripe-II
state has the lowest ground state energy with an extrapolated value of
$E_{\infty}=-2.8500\pm 0.0026$ in the thermodynamic limit.

This strip-II state breaks the lattice $C_{6}$ rotational symmetry to a
$C_{2}$ rotational symmetry and has $2\times{}1$ unit cells. Since the SU(4)
spin-orbital rotational symmetry is respected, such a state is indeed a
nematic quantum spin-orbital liquid state. The low energy excitations on top
of the ground state are gapless and can be represented by a parton Fermi
surface consisting of open orbits in the Brillouin zone. All of these results
are in excellent agreement with previous results obtained from the Gutzwiller-
boosted DMRG [53].

###### Acknowledgements.

We would like to thank Rong-Yang Sun and Hong-Hao Tu for their helpful
discussions. This work is supported in part by National Key Research and
Development Program of China (No. 2022YFA1403403), National Natural Science
Foundation of China (No.12274441, 12034004), and the European Research Council
(ERC) under the European Unions Horizon 2020 research and innovation program
(Grant agreement No. 771537).

##  Appendix A Stochastic reconfiguration method in VMC calculations

To perform the variational Monte Carlo (VMC) optimization, a stochastic
reconfiguration (SR) method has been adopted. It simulates the effect of the
Hessian matrix with a more cost-effective approach, which we now review
briefly.

When carrying out the VMC optimization, one key step is the computation of the
expectation value of Hamiltonian and its gradient with respect to the
variational parameters. Without loss of generality, the variational state
$|\Psi\rangle$ can be constructed by the wave function $|\Psi(C)\rangle$ in an
orthonormal basis $|C\rangle$, written as

$|\Psi\rangle=\sum_{C}\Psi(C)|C\rangle.$ (A1)

Then, the variational energy is given by

$E=\frac{\langle\Psi|H|\Psi\rangle}{\langle\Psi\mid\Psi\rangle}=\frac{\sum_{C}|\Psi(C)|^{2}E_{loc}(C)}{\sum_{C}|\Psi(C)|^{2}}.$
(A2)

The local energy $E_{loc}$ is sampled and its average value is calculated by:

$E_{loc}(C)=\sum_{C^{\prime}}\frac{\langle
C|H|C^{\prime}\rangle\Psi(C^{\prime})}{\Psi(C)}.$ (A3)

To search for the minimum energy of $E(t_{1},\cdots,t_{12})$ in which a value
taking condition $\sum_{n=1}^{12}t_{n}=12$ has been imposed. For convenience,
the series of parameters in wave function
$|\Psi(t_{1},t_{2},\cdots,t_{12})\rangle$ have been abbreviated as
$\\{\theta\\}$. The gradient of the variational energy with respect to the
variational parameters $\\{\theta\\}$ is given by

$\begin{split}\nabla_{\theta}E(\theta)=&2\langle\nabla_{\theta}(ln\Psi(C))E_{loc}(C)\rangle\\\
&-2\langle
E_{loc}(C)\rangle\langle\nabla_{\theta}(ln\Psi(C))\rangle,\end{split}$ (A4)

which is obtained by the derivative with respect to $\\{\theta\\}$ of the
following formula

$E(\theta)=\frac{\langle\Psi(\theta)|H|\Psi(\theta)\rangle}{\langle\Psi(\theta)|\Psi(\theta)\rangle}.$
(A5)

In the formula (A4), the two quantities can be calculated by standard Monte
Carlo sampling on the distribution generated by $|\Psi(\theta)|^{2}$, and the
quantity $\langle\nabla_{\theta}(ln\Psi(C))\rangle$ will play a key role in
the calculation of the gradient.

To obtain the adjustment direction and the step length of the $\\{\theta\\}$,
an approximate Hessian matrix
$\langle\Psi(\theta)|H|\dot{\Psi}(\theta)\rangle$ requires the construction

$\frac{\partial^{2}}{\partial\theta\partial\theta^{\prime}}\frac{\langle\Psi(\theta)|\Psi(\theta^{\prime})\rangle}{\langle\Psi(\theta)|\Psi(\theta)\rangle}.$
(A6)

As the essence of the SR method, the second derivative of energy with respect
to parameters $\\{\theta\\}$ can be obtained by matrix (A6). For imitating the
effect of the Hessian matrix, a positive-definite Hermitian matrix $S$
generated from the metric of the variational state in the variational space is
proposed. More specifically, the matrix element is given by

$\begin{split}S_{ij}=&\langle\nabla_{\theta}[ln\Psi(C)]_{i}\nabla_{\theta}[ln\Psi(C)]_{j}\rangle-\langle\nabla_{\theta}[ln\Psi(C)]_{i}\rangle\langle\nabla_{\theta}[ln\Psi(C)]_{j}\rangle,\end{split}$
(A7)

where the index $i$ is with respect to the $i$-th variational paramter.

Thanks to the quantity $\langle\nabla_{\theta}(ln\Psi(C))\rangle$ derived from
Eq.(A4) above, an approximate Hessian matrix $S$ can be obtained from Eq.
(A7), Therefore, the direction and step length for performing parameters
optimization can be found as follows,

$\Delta\vec{\theta}=-\frac{1}{2\kappa}S^{-1}\vec{g}.$ (A8)

Here $\kappa$ is an empirical positive parameter, which is adjusted manually
and defines the step length of parameter optimization.

##  Appendix B VMC calculations for four Gutzwiller projected states listed in
Table 1

Table A1: The optimized ground states for four types of Gutzwiller projected states listed in Table 1. Here $E_{L}$ is the optimized ground state energy per site on an XC($L\times{}L$) torus with $L=12,16,24$. The variational parameters $\delta_{1}$, $\delta_{2}$ and $\epsilon$ are defined in Eqs. (10), (11), (12), (13), and (14), respectively.  | $L=12$ | $L=16$ | $L=24$
---|---|---|---
| $E_{L}$ | $\delta_{1},\delta_{2},\epsilon$ | $E_{L}$ | $\delta_{1},\delta_{2},\epsilon$ | $E_{L}$ | $\delta_{1},\delta_{2},\epsilon$
(a) stripe-I state | $-2.8091\pm 0.0014$ |  $\delta_{1}\approx 0.174$ | $-2.8276\pm 0.0021$ |  $\delta_{1}\approx 0.156$ | $-2.8377\pm 0.0024$ |  $\delta_{1}\approx 0.145$
(b) stripe-II state | $-2.8174\pm 0.0015$ |  $\delta_{1}\approx 0.181$ $\epsilon\approx{}0$ | $-2.8369\pm 0.0021$ |  $\delta_{1}\approx 0.161$ $\epsilon\approx{}0$ | $-2.8462\pm 0.0024$ |  $\delta_{1}\approx 0.151$ $\epsilon\approx{}0$
(c) stripe-III state | $-2.8015\pm 0.0014$ |  $\delta_{1}\approx 0.177$ | $-2.8210\pm 0.0020$ |  $\delta_{1}\approx 0.156$ | $-2.8312\pm 0.0025$ |  $\delta_{1}\approx 0.145$
(d) plaquette state | $-2.7401\pm 0.0010$ |  $\delta_{1}\approx 0.218$ $\delta_{2}\approx 0.217$ | $-2.7600\pm 0.0018$ |  $\delta_{1}\approx\delta_{2}\approx 0.203$ | $-2.7670\pm 0.0020$ |  $\delta_{1}\approx\delta_{2}\approx 0.208$

In the main text, using the SR method, we have searched for local energy
minima on XC($L\times{L}$) tori up to $L=24$ (see Fig. 3 and Table 2) with the
generic 12-parameter $2\times{}2$ VBC states given in Eq. (10). In this
appendix, we perform standard VMC calculations and search for local energy
minima for four Gutzwiller projected states on XC($L\times{}L$) tori with
$L=12,16,24$. These Gutzwiller projected states are characterized by the one
or two variational parameters, as listed in Table 1, namely, $\delta_{1}$ in
Eq. (12) and Eq. (14) [stripe-I and stripe-III states], $\delta_{1}$ and
$\epsilon$ in Eq. (13) [stripe-II state], and $\delta_{1}$ and $\delta_{2}$ in
Eq. (15) [plaquette state].

The optimized ground state energy and corresponding variational parameters are
given in Table A1, which are consistent with the 12-parameter optimization
results shown in Table 2 and Fig. 3 in the main text.

## References

  * Anderson [1973] P. Anderson, Resonating valence bonds: A new kind of insulator?, Mater. Res. Bull. 8, 153 (1973).
  * Lee [2008] P. A. Lee, An end to the drought of quantum spin liquids, Science 321, 1306 (2008).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2016).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A field guide to spin liquids, Annu. Rev. Condens. Matter Phys. 10, 451 (2019).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367 (2020).
  * Huse and Elser [1988] D. A. Huse and V. Elser, Simple variational wave functions for two-dimensional heisenberg spin-½ antiferromagnets, Phys. Rev. Lett. 60, 2531 (1988).
  * Jolicoeur and Le Guillou [1989] T. Jolicoeur and J. C. Le Guillou, Spin-wave results for the triangular heisenberg antiferromagnet, Phys. Rev. B 40, 2727 (1989).
  * White and Chernyshev [2007] S. R. White and A. L. Chernyshev, Neél order in square and triangular lattice heisenberg models, Phys. Rev. Lett. 99, 127004 (2007).
  * Jin and Zhou [2022] H.-K. Jin and Y. Zhou, Material search for quantum spin liquids on the simplest frustrated lattice, Chin. Phys. Lett. 39, 050102 (2022).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Spin-liquid ground state of the s= 1/2 kagome heisenberg antiferromagnet, Science 332, 1173 (2011).
  * Depenbrock  _et al._ [2012] S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Nature of the spin-liquid ground state of the $s=1/2$ heisenberg model on the kagome lattice, Phys. Rev. Lett. 109, 067201 (2012).
  * Liao  _et al._ [2017] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, Gapless spin-liquid ground state in the $s=1/2$ kagome antiferromagnet, Phys. Rev. Lett. 118, 137202 (2017).
  * He  _et al._ [2017] Y.-C. He, M. P. Zaletel, M. Oshikawa, and F. Pollmann, Signatures of dirac cones in a dmrg study of the kagome heisenberg model, Phys. Rev. X 7, 031020 (2017).
  * Block  _et al._ [2011] M. S. Block, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Spin bose-metal and valence bond solid phases in a spin-$1/2$ model with ring exchanges on a four-leg triangular ladder, Phys. Rev. Lett. 106, 157202 (2011).
  * Mishmash  _et al._ [2013] R. V. Mishmash, J. R. Garrison, S. Bieri, and C. Xu, Theory of a competitive spin liquid state for weak mott insulators on the triangular lattice, Phys. Rev. Lett. 111, 157203 (2013).
  * Zhu and White [2015] Z. Zhu and S. R. White, Spin liquid phase of the $s=\frac{1}{2}\phantom{\rule{4.0pt}{0.0pt}}{J}_{1}-{J}_{2}$ heisenberg model on the triangular lattice, Phys. Rev. B 92, 041105 (2015).
  * Hu  _et al._ [2015] W.-J. Hu, S.-S. Gong, W. Zhu, and D. N. Sheng, Competing spin-liquid states in the spin-$\frac{1}{2}$ heisenberg model on the triangular lattice, Phys. Rev. B 92, 140403 (2015).
  * Wietek and Läuchli [2017] A. Wietek and A. M. Läuchli, Chiral spin liquid and quantum criticality in extended $s=\frac{1}{2}$ heisenberg models on the triangular lattice, Phys. Rev. B 95, 035141 (2017).
  * Gong  _et al._ [2017] S.-S. Gong, W. Zhu, J.-X. Zhu, D. N. Sheng, and K. Yang, Global phase diagram and quantum spin liquids in a spin-$\frac{1}{2}$ triangular antiferromagnet, Phys. Rev. B 96, 075116 (2017).
  * Saadatmand and McCulloch [2017] S. N. Saadatmand and I. P. McCulloch, Detection and characterization of symmetry-broken long-range orders in the spin-$\frac{1}{2}$ triangular heisenberg model, Phys. Rev. B 96, 075117 (2017).
  * Hu  _et al._ [2019] S. Hu, W. Zhu, S. Eggert, and Y.-C. He, Dirac spin liquid on the spin-$1/2$ triangular heisenberg antiferromagnet, Phys. Rev. Lett. 123, 207203 (2019).
  * Feiner  _et al._ [1997] L. F. Feiner, A. M. Oleś, and J. Zaanen, Quantum melting of magnetic order due to orbital fluctuations, Phys. Rev. Lett. 78, 2799 (1997).
  * Kugel and Khomskii [1982] K. I. Kugel and D. I. Khomskii, The jahn-teller effect and magnetism: transition metal compounds, Soviet Physics Uspekhi 25, 231 (1982).
  * Li  _et al._ [1998] Y. Q. Li, M. Ma, D. N. Shi, and F. C. Zhang, Su(4) theory for spin systems with orbital degeneracy, Phys. Rev. Lett. 81, 3527 (1998).
  * Pati  _et al._ [1998] S. K. Pati, R. R. P. Singh, and D. I. Khomskii, Alternating spin and orbital dimerization and spin-gap formation in coupled spin-orbital systems, Phys. Rev. Lett. 81, 5406 (1998).
  * Khaliullin and Maekawa [2000] G. Khaliullin and S. Maekawa, Orbital liquid in three-dimensional mott insulator: ${\mathrm{latio}}_{3}$, Phys. Rev. Lett. 85, 3950 (2000).
  * Tokura and Nagaosa [2000] Y. Tokura and N. Nagaosa, Orbital physics in transition-metal oxides, Science 288, 462 (2000).
  * Reynaud  _et al._ [2001] F. Reynaud, D. Mertz, F. Celestini, J.-M. Debierre, A. M. Ghorayeb, P. Simon, A. Stepanov, J. Voiron, and C. Delmas, Orbital frustration at the origin of the magnetic behavior in ${\mathrm{linio}}_{2}$, Phys. Rev. Lett. 86, 3638 (2001).
  * Yamada  _et al._ [2018] M. G. Yamada, M. Oshikawa, and G. Jackeli, Emergent $\mathrm{SU}(4)$ Symmetry in $\alpha\text{$-$}{\mathrm{ZrCl}}_{3}$ and Crystalline Spin-Orbital Liquids, Phys. Rev. Lett. 121, 097201 (2018).
  * Yamada  _et al._ [2021] M. G. Yamada, M. Oshikawa, and G. Jackeli, Su(4)-symmetric quantum spin-orbital liquids on various lattices, Phys. Rev. B 104, 224436 (2021).
  * Sutherland [1975] B. Sutherland, Model for a multicomponent quantum system, Phys. Rev. B 12, 3795 (1975).
  * Li  _et al._ [1999] Y.-Q. Li, M. Ma, D.-N. Shi, and F.-C. Zhang, Ground state and excitations of a spin chain with orbital degeneracy, Phys. Rev. B 60, 12781 (1999).
  * Affleck [1986] I. Affleck, Exact critical exponents for quantum spin chains, non-linear $\sigma$-models at $\theta=\pi$ and the quantum hall effect, Nucl. Phys. B 265, 409 (1986).
  * Azaria  _et al._ [1999] P. Azaria, A. O. Gogolin, P. Lecheminant, and A. A. Nersesyan, One-dimensional su(4) spin-orbital model: A low-energy effective theory, Phys. Rev. Lett. 83, 624 (1999).
  * Führinger  _et al._ [2008] M. Führinger, S. Rachel, R. Thomale, M. Greiter, and P. Schmitteckert, Dmrg studies of critical su (n) spin chains, Annalen der Physik 17, 922 (2008).
  * van den Bossche  _et al._ [2001] M. van den Bossche, P. Azaria, P. Lecheminant, and F. Mila, Spontaneous plaquette formation in the su(4) spin-orbital ladder, Phys. Rev. Lett. 86, 4124 (2001).
  * Weichselbaum  _et al._ [2018] A. Weichselbaum, S. Capponi, P. Lecheminant, A. M. Tsvelik, and A. M. Läuchli, Unified phase diagram of antiferromagnetic su($n$) spin ladders, Phys. Rev. B 98, 085104 (2018).
  * Chen  _et al._ [2005] S. Chen, C. Wu, S.-C. Zhang, and Y. Wang, Exact spontaneous plaquette ground states for high-spin ladder models, Phys. Rev. B 72, 214428 (2005).
  * Corboz  _et al._ [2012] P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Spin-orbital quantum liquid on the honeycomb lattice, Phys. Rev. X 2, 041013 (2012).
  * Natori  _et al._ [2019] W. M. H. Natori, R. Nutakki, R. G. Pereira, and E. C. Andrade, Su(4) heisenberg model on the honeycomb lattice with exchange-frustrated perturbations: Implications for twistronics and mott insulators, Phys. Rev. B 100, 205131 (2019).
  * Jin  _et al._ [2023] H.-K. Jin, W. M. H. Natori, and J. Knolle, Twisting the dirac cones of the su(4) spin-orbital liquid on the honeycomb lattice, Phys. Rev. B 107, L180401 (2023).
  * van den Bossche  _et al._ [2000] M. van den Bossche, F. C. Zhang, and F. Mila, Plaquette ground state in the two-dimensional su (4) spin-orbital model, Eur. Phys. J. B 17, 367 (2000).
  * Wang and Vishwanath [2009] F. Wang and A. Vishwanath, ${\text{z}}_{2}$ spin-orbital liquid state in the square lattice kugel-khomskii model, Phys. Rev. B 80, 064413 (2009).
  * Corboz  _et al._ [2011] P. Corboz, A. M. Läuchli, K. Penc, M. Troyer, and F. Mila, Simultaneous dimerization and su(4) symmetry breaking of 4-color fermions on the square lattice, Phys. Rev. Lett. 107, 215301 (2011).
  * Penc  _et al._ [2003] K. Penc, M. Mambrini, P. Fazekas, and F. Mila, Quantum phase transition in the su(4) spin-orbital model on the triangular lattice, Phys. Rev. B 68, 012408 (2003).
  * Keselman  _et al._ [2020] A. Keselman, B. Bauer, C. Xu, and C.-M. Jian, Emergent fermi surface in a triangular-lattice su(4) quantum antiferromagnet, Phys. Rev. Lett. 125, 117202 (2020).
  * Zhang  _et al._ [2021] Y.-H. Zhang, D. N. Sheng, and A. Vishwanath, Su(4) chiral spin liquid, exciton supersolid, and electric detection in moiré bilayers, Phys. Rev. Lett. 127, 247701 (2021).
  * Yao  _et al._ [2021] X.-P. Yao, Y. Gao, and G. Chen, Topological chiral spin liquids and competing states in triangular lattice $\mathrm{SU}(n)$ mott insulators, Phys. Rev. Research 3, 023138 (2021).
  * Jin  _et al._ [2021] H.-K. Jin, H.-H. Tu, and Y. Zhou, Density matrix renormalization group boosted by gutzwiller projected wave functions, Phys. Rev. B 104, L020409 (2021).
  * Jin  _et al._ [2022] H.-K. Jin, R.-Y. Sun, H.-H. Tu, and Y. Zhou, Unveiling a critical stripy state in the triangular-lattice su (4) spin-orbital model, Science Bulletin 67, 918 (2022).
  * Georgi [1999] H. Georgi,  _Lie algebras in particle physics_ (Perseus Books, 1999).
  * Sorella [1998] S. Sorella, Green function monte carlo with stochastic reconfiguration, Phys. Rev. Lett. 80, 4558 (1998).
  * Sorella and Capriotti [2000] S. Sorella and L. Capriotti, Green function monte carlo with stochastic reconfiguration: An effective remedy for the sign problem, Phys. Rev. B 61, 2599 (2000).
  * Manousakis [1991] E. Manousakis, The spin-½ heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides, Rev. Mod. Phys. 63, 1 (1991).
  * Gross  _et al._ [1989] M. Gross, E. Sánchez-Velasco, and E. Siggia, Ground-state properties of the two-dimensional antiferromagnetic heisenberg model, Phys. Rev. B 39, 2484 (1989).
  * Neuberger and Ziman [1989] H. Neuberger and T. Ziman, Finite-size effects in heisenberg antiferromagnets, Phys. Rev. B 39, 2608 (1989).


Progressive Semantic Segmentation

Chuong Huynh^1   Anh Tuan Tran^1,2   Khoa Luu^1,3   Minh Hoai^1,4

^1VinAI Research, Hanoi, Vietnam,  ^2VinUniversity, Hanoi, Vietnam

^3University of Arkansas, Fayetteville, AR 72701, USA

^4Stony Brook University, Stony Brook, NY 11790, USA

The objective of this work is to segment high-resolution images without overloading GPU memory usage or losing the fine details in the output segmentation map. The memory constraint means that we must either downsample the big image or divide the image into local patches for separate processing. However, the former approach would lose the fine details, while the latter can be ambiguous due to the lack of a global picture.
In this work, we present MagNet, a multi-scale framework that resolves local ambiguity by looking at the image at multiple magnification levels. MagNet has multiple processing stages, where each stage corresponds to a magnification level, and the output of one stage is fed into the next stage for coarse-to-fine information propagation. Each stage analyzes the image at a higher resolution than the previous stage, recovering the previously lost details due to the lossy downsampling step, and the segmentation output is progressively refined through the processing stages. Experiments on three high-resolution datasets of urban views, aerial scenes, and medical images show that MagNet consistently outperforms the state-of-the-art methods by a significant margin. Code is available at <https://github.com/VinAIResearch/MagNet>.

§ INTRODUCTION

[](a) Input image

[](b) Ground truth

[](c) Downsampling

[](d) Patch processing

[](e) GLNet

[](f) DenseCRF

[](g) PointRend

[](h) MagNet (Proposed)

Comparing several semantic segmentation and refinement approaches on a high-resolution input image. Downsampling loses fine details, while Patch Processing wrongly classifies local patches due to the lack of the global context. The collaborative global-local network GLNet fails due to the large discrepancy between the global and local branches. Post-processing and refinement methods such as DenseCRF and PointRend  can only correct small mistakes due to local inconsistency. MagNet outperforms other methods, thanks to a novel multi-scale segmentation and refinement framework.
Best viewed in color.

The current state-of-the-art (SOTA) semantic image segmentation techniques <cit.> are based on deep learning, where a convolutional neural network (CNN) takes an input image and outputs a segmentation map. Most of the existing techniques, however, assume that the entire segmentation process can be performed with a single feed-forward pass of the input image and the entire process can be fitted into GPU memory. Unfortunately, most existing techniques cannot handle high-resolution input images due to memory and other computational constraints.
One approach to handle a large input image is to downsample it, but this results in a low-resolution segmentation map, which is not adequate for applications that require high-resolution output with fine details <cit.>, e.g., for tracking the progression of malignant lesion <cit.>. Another approach to handle a large input image is to divide the image into small patches and process each patch independently. This approach, however, does not take into account the global information <cit.> that is needed to resolve ambiguity in local patches. The limitations of these two approaches are illustrated in fig:teaser(c) & (d).

One way to address the limitations of the two aforementioned approaches is to combine them, i.e., to fuse global and local segmentation processes. On the one hand, the global view of the entire image can be used to resolve the ambiguity in the appearance of local patches. On the other hand, by analyzing local patches, we can refine the segmentation boundaries and recover the lost details due to the downsampling procedure of the global segmentation process. This approach has been successfully demonstrated recently by the Global-Local Network (GLNet) <cit.>. However, given an input image with ultra-high resolution, there is a huge gap between the scale of the whole image and the scale of the local patches. This will lead to contrasting output segmentation maps, and it will be difficult to combine and reconcile differences with a single feed-forward processing stage (see fig:teasere); the difficulty of this combination task is analogous to constructing a single-span bridge across a wide river.

To bridge the gap between the two extreme ends of the scale space, we propose to consider multiple scales in between. We introduce a novel multi-scale framework where the output segmentation map will be progressively refined as the image is analyzed from the coarsest to the finest scale.
The core of our framework is a refinement module that can use one segmentation map to refine another. This refinement module is used at every stage of our multi-scale processing pipeline to refine the output segmentation map at its most uncertain locations. Our framework can integrate global contextual cues to produce more accurate segmentation, and it can output high-resolution detailed segmentation maps under a memory constraint. fig:teaser shows the result of MagNet and compares it with other segmentation methods, including the recently proposed PointRend <cit.> method that seeks to refine only at the most uncertain pixels.

§ RELATED WORK

to @X[5.2,l] X[1.9,c] X[1,c] X[1,c] X[1,c] X[1.2,c] X[1.6,c] X[1.9,c] X[1.2,c] X[1.7,c] X[1,c]@

CRF <cit.>
PSP <cit.>

Deep learning

Using high-resolution input

Multi-scale processing

Can recognize new objects

Partly refinement

Summary of key features of various semantic segmentation refinement approaches.

Multi-scale, multi-stage, context aggregation. The combination of multiple scale levels helps the network aggregate different fields of view and provides more context to each pixel <cit.>. ICNet <cit.> used a cascaded architecture for feature maps of different downsampled inputs, while RefineNet <cit.> fused upsampled outputs of the branches that handled different low-resolution inputs. Feature Pyramid Network (FPN) <cit.> upsampled feature maps in different scales and aggregated them with the output of low layers. The dilated convolution and Atrous Spatial Pyramid Pooling (ASPP) module in DeepLab <cit.> enlarged the receptive field, created a connection between far-apart pixels. The same effect was achieved by PSPNet <cit.>, which combined different scaling feature maps to enlarge the receptive fields. High-resolution Net (HRNet) <cit.> proposed another scale fusion schema where a new branch with a larger receptive field was added after each stage. Attention technique was also used in many recent approaches <cit.> to add more global information to every single point.

Another approach to handle high-resolution images is to use multi-stage networks, where images are segmented after several stages or sub-networks. <cit.> proposed Hierarchical Auto-Zoom Net, a strategy to scale the field of view when sliding the view through a big image. For ultra-high resolution images, <cit.> solved the imbalance between background and foreground by predicting whether the whole patch contained foreground pixels or not before segmenting it.

Propagating global information to local patches is a promising approach to deal with high-resolution images. ParseNet <cit.> pooled global context to local field of view to have more information. BiSeNet <cit.> included one more branch for global pooling and the global context would be added to the feature map at the last stage. Although those methods were efficient, they consumed a huge amount of GPU memory for ultra-high resolution images. <cit.> proposed a method for super-high-resolution images, using independent multi-scale networks and an adaptive weight generator. The outputs of network members were combined with corresponding trained weights to produce the final output, but there was no knowledge sharing between network branches. Unlike <cit.>, <cit.> contained two sub-networks with shared information, where the global branch took the downsampled images to extract global context and the local branch took patches and corresponding global features to improve the details of high-resolution images. However, due to the ad-hoc combination of the global and local branches, it was difficult to extend to more than two scales.  Moreover, in our experience, training the local network was difficult due to the domination of the strong global branch.

Segmentation refinement. There were several approaches to improve the segmentation outputs with post-processing. One approach was to use classical methods such as Conditional Random Fields (CRFs) <cit.> or Guided Filter (GF) <cit.> on the segmentation mask produced by deep learning networks. However, these methods were slow and the improvement was incremental. The inference speed could be improved with a deep learning version of Guided Filter (DGF) <cit.>. Another approach for post-processing was to use deep networks. Iterative Instance Segmentation (ISS) <cit.> refined the output by repeatedly passing the input image and the segmentation map through a refinement module several times. This method was based on self-reflection, the input image to each refinement stage was the same. Like ISS, CascadedPSP <cit.> used the same refinement scheme but the resolution of the input at each refinement stage was different. However, the wrong prediction at any middle stage could significantly affect the performance of later steps. Some methods aimed to refine parts of the output only, such as pixels at the boundaries <cit.> or pixels at uncertain locations (PointRend <cit.>). However, boundary refinement methods <cit.> failed to recover tiny objects, while PointRend <cit.> only used the local context for refinement.
Furthermore, because the input of the PointRend was the high-level features of a deep network, it must be trained specifically for each segmentation backbone. In this paper, we propose a modular framework for having any number of scale levels. It is simple but effective for refining a coarse segmentation output, being able to keep the overall structure of the coarse segmentation output while adding more details after each stage without suffering the domination problems. tab:method_sum compares the key features of different methods.

§ MAGNET

Overview of our proposed MagNet. The segmentation network produces the scale-specific prediction while the refinement module selectively refines the coarse prediction from previous stages based on that local prediction.

We now describe MagNet, a multi-scale segmentation framework for ultra-high resolution images. It is a multi-stage network architecture, where each stage of the network corresponds to a particular scale. An input image will be inspected at multiple scales, from the coarsest to the finest.

The core of our framework is a segmentation module and a refinement module, which are used at every processing stage. At each stage, the inputs to the refinement module are two segmentation maps: (1) the cumulative result from the previous stages, and (2) the result obtained by running the segmentation module at and only at the current scale. The objective of the refinement module is to use the latter segmentation map to refine the former one, at selective locations determined based on the uncertainty of two estimated segmentation maps.

In our framework, the segmentation module can be any segmentation backbone, as long as it can output a segmentation map with uncertainty estimates. The refinement module is agnostic to the segmentation backbone; it can be trained with one backbone and used with another. In the following, we will describe the multi-stage processing pipeline and the refinement module in details.

§.§ Multi-stage processing pipeline

The architecture and processing pipeline of MagNet is depicted in fig:overall_structure. There is one segmentation module and one refinement module, which are used repeatedly in $m$ processing stages, where $m$ is a hyper-parameter for the number of scales that we want to analyze. We use $s$ to denote the processing stage, where $s=1$ corresponds to the coarsest scale and $s=m$ corresponds to the finest scale. Let $X \in \mathbb{R}^{H{\times}W{\times}3}$ be an input image, where $H, W$ are the height and width of the image. We consider the case when $H$ and $W$ are too big for image $X$ to be processed without downsampling, and let $h$ and $w$ be the largest  (or desirable) sizes that can be handled by the segmentation module. We use $h^s$ and $w^s$ to denote the height and width for the scale level $s$. We determine the scale levels so that they span the entire scale space: $H = h^1 > \cdots > h^m = h$ and $W = w^1 > \cdots > w^m = w$.

For a particular scale level $s$, we divide the input image $X$ into patches of size $\htimesw{h^s}{w^s}$ and perform semantic segmentation on these patches. The locations of these patches are defined by a set of rectangular windows, and let $\mP^s$ denote the set of these windows: $\mP^s = \{\p| \p = (x, y, w^s, h^s) \}$, where each window is specified by the top-left corner, width, and height. As the scale level $s$ increases, the width and height of the rectangular windows decrease, but the cardinality of $\mP^s$ increases. For a particular window $\p$, we will use $X_{\p}$ to denote the image patch extracted at the window $\p$.

Our network will take an image $X \in \mathbb{R}^{H{\times}W{\times}3}$ and produce a sequence of segmentation maps $Y^1, \cdots, Y^m \in \mathbb{R}^{H{\times}W{\times}C}$, where $C$ is the number of semantic categories in consideration. At stage $s$, we first determine the set of rectangular windows  $\mP^s$ for patch division and refine the segmentation map of each patch. Specifically, for each window $\p \in \mP^s$, do:

* Extract the image patch $X^s_{\p}$ and previous segmentation output $Y^{s-1}_\p$ defined by the window $\p$. The height and width of these tensors are $h^s$ and $w^s$.

* Downsample $X^s_{\p}$ and $Y^{s-1}_\p$ so that the new height and width are $h$ and $w$, which are the size that can be fitted into GPU memory and be processed by the segmentation and refinement modules. Let $\bar{X}^s_{\p}$ and $\bar{Y}^{s-1}_{\p}$ denote the downsampled tensors.

* With $\bar{X}^s_{\p}$ as the input, use the segmentation module to obtain the scale-specific segmentation map $\bar{O}^s_\p$.

* With $\bar{Y}^{s-1}_{\p}$ and $\bar{O}^s_\p$ as the input, use the refinement module to refine $\bar{Y}^{s-1}_{\p}$ to obtain $\bar{Y}^s_{\p}$ (See sec:refinement_module).

* Upsample $\bar{Y}^s_{\p}$ to get $Y^s_{\p}$ of size $h^s{\times}w^s{\times}C$.

These processing steps are illustrated in fig:overall_structure.

§.§ Refinement module

The refinement module is a core component of our framework, which is used to refine the individual patches of a segmentation map at every processing stage of our pipeline. The input to this module is two segmentation maps of size $h{\times}w{\times}C$: (1) the scale-cumulative segmentation map $Y$, from all previous scales, and (2) the scale-specific segmentation map $O$, from the current scale. The output of the module is the updated scale-cumulative segmentation map.

The overview of the refinement module. The cumulative segmentation $Y$ is partly replaced with the scale-specific segmentation map $O$ based on the score $\mQ$.

fig:refinement_module depicts the refinement process, which contains the following steps. First, using a small network with $Y$ and $O$ as the input, we obtain an initial combined segmentation map $R$. We then calculate the prediction uncertainty maps for both $Y$ and $R$.  Specifically, for each pixel of $Y$, the prediction confidence at this location is defined as the absolute difference between the highest probability value and the second-highest value (among the $C$ probability values for $C$ classes). The uncertainty score is then computed based on the confidence score such that the two scores must add up to one. Similarly, we can compute the prediction uncertainty map for $R$. Let $Y^u$ and $R^u$ denote the prediction uncertainty maps for $Y$ and $R$ respectively.

Next, we will describe how to use two uncertainty maps to select $k$ locations of $Y$ for refinement. While only one uncertainty prediction is used in the previous work PointRend <cit.>, we extend this approach to have a better selection strategy. These are the locations where $Y$ is uncertain about its prediction, while $R$ is certain about its prediction. The score map for ranking the pixels is calculated as $\mQ = \mF(Y^u \odot (1 - R^u))$, where $\odot$ denotes point-wise multiplication and $\mF$ denotes median blurring to smooth the score map. Empirical comparison between the effect of each element in the formula can be found in sec:ablation_selection.

§.§ MagNet-Fast

There is trade-off between the accuracy and the run-time efficiency of our framework. One way to reduce the running time is to decrease the number of scales to process. Another approach is to perform segmentation and refinement on a subset of image patches at each scale level. MagNet-Fast combines these two approaches when it runs on a smaller number of scales and only selects the patches with the highest prediction uncertainty $Y^u$ for refinement. In MagNet-Fast, the total number of image patches that need to be fed into the segmentation module might be even smaller than the number of image patches used in the patch processing approach. Moreover, MagNet-Fast can leverage both global context and detailed information for segmentation, leading to superior results as will be seen in our experiments.

§ EXPERIMENTS

We evaluated the performance of MagNet on three high resolution datasets: Cityscapes <cit.>, DeepGlobe <cit.>, and Gleason <cit.>. Some information about these datasets is listed in tab:datasets. The number of pixels of each image is from 2 to 25 million. In this section, we present experiments comparing MagNet with other state-of-the-art frameworks in semantic segmentation and also describe some ablation studies on Cityscapes.

to @X[3,l]X[3,l]X[2,c]X[2,c]@
Dataset     Content     Resolution     No.classes

Cityscapes <cit.>     urban scene     $2048{\times}1024$     19

DeepGlobe <cit.>     aerial scene     $2448{\times}2448$     6

Gleason <cit.>     histopathology     $5000{\times}5000$     4

Details of high-resolution datasets used to evaluate our framework. All images have from 2 to 25 million pixels with a lot of details.

§.§ Implementation details

The two residual blocks are trained to refine the segmentation at each scale. This module outputs the same size $h{\times}w$ as the input.

Architecture of the refinement module. fig:refinement_network depicts the architecture of the refinement module used in all experiments. The main components are the two residual blocks. With the input of size $h{\times}w{\times}2C$, the refinement module produces the output of size $h{\times}w{\times}C$.

Training. For each dataset, we trained a state-of-the-art segmentation model on the downsampled images and a refinement module to refine the coarse output on sliced images. While training the refinement module, we randomly extracted image patches and also applied the following data augmentation: rotation, and horizontal and vertical flipping. We used SGD optimizer with momentum $0.9$, decayed weight $5{\times}10^{-4}$, and initial learning rate $10^{-3}$. We trained the refinement module for 50 epochs, and the learning rate was decreased tenfold at epoch 20, 30, 40, and 45. We used cross-entropy as the loss function for training segmentation and refinement modules. We implemented MagNet using PyTorch <cit.> starting from the public implementation of HRNet-OCR <cit.>. We use a batch size of 16 for training on a DGX-1 workstation with Tesla V100 GPUs.

Testing. During inference, at each scale, we extracted non-overlapping patches for processing. The evaluation metric is mean Intersection over Union (mIoU). For memory and speed comparison, we ran benchmarking experiments on a machine with an Intel i7 CPU and an RTX2080Ti GPU.

§.§ Experiments on the Cityscapes dataset

Cityscapes is a dataset of high-resolution urban scenes, containing images of size $2048{\times}1024$ pixels. The task is to segment objects in videos captured by auto cameras. There are two kinds of data, with coarse and fine labels. In our experiments, we used the fine-labeled dataset, with 2,975 training images and 500 testing images.

We considered four possible scale levels ($m{=}4$), corresponding to patch sizes of $2048{\times}1024$, $1024{\times}512$, $512{\times}256$, and $256{\times}128$. The size of the input to the segmentation module was always $256{\times}128$ to satisfy the memory constraint, so any larger patch would be downsampled.

§.§.§ Benefits of multiple scale levels

to @X[3,l]X[1,c]X[1,c]@
Refinement steps     mIoU(%)     Time(s)

256      63.23     0.03

256$\rightarrow$512     65.73     0.19

256$\rightarrow$1024     65.23     0.61

256$\rightarrow$2048     65.21     2.24

256$\rightarrow$512$\rightarrow$2048     67.13     2.38

256$\rightarrow$1024$\rightarrow$2048     66.95     2.79

256$\rightarrow$512$\rightarrow$1024$\rightarrow$2048     67.57     2.93

Performance of MagNet on Cityscapes with and without intermediate scale levels. It is essential to have the intermediate scales.

The visualization of segmentation output through each processing stage on Cityscapes dataset. The first one is the image. Others in the first row are the selective points to be refined (red color). The second row is the segmentation output. The third row is the errors comparing with the ground-truth. Best viewed in color.

tab:cityscapes_scales shows the results of MagNet for a different number of scales.  While the direct refinement from the lowest to highest scale improves about 2%  mIoU, from 63.23% to 65.23%, adding the two intermediate scales between the smallest and largest scales improve the performance by 4.34% mIoU. Qualitative improvements through each processing stage can be observed in fig:step_refinement. After each stage, the errors decrease and the segmentation masks become finer.

§.§.§ Comparing segmentation approaches

The IoU cumulative distribution

Visualization of refinement methods

Our methods outperform other refinement frameworks on the Cityscapes dataset. (a) The cumulative distribution of mIoU of each image on the dataset (lower is better). The MagNet and MagNet-Fast achieve the best result among others. (b) Some segmentation results of refinement frameworks and our MagNet. The mIoU numbers are below the images. More tiny objects are recognized and boundaries are refined better with MagNet and MagNet-Fast. Best viewed in a digital device with magnification.

to @X[3.5,l]X[3,l]X[2.5,l]X[1.5,r]@
Model     mIoU(%)     Mem.(MB)     Time(s)

Patching     52.19     1575     1.77

Downsampling      63.23     1575     0.02

DenseCRF <cit.>     62.95 0.28     1575     26.02

DGF <cit.>     63.33 0.10     1727 152     0.32

PointRend <cit.>     64.39 1.16     2033 458     0.14

SegFix <cit.>     65.83 2.60     2961 1386     0.38

MagNet-Fast     66.91 3.68     2007 432     0.32

MagNet     67.57 4.34     2007 432     2.93

Performance of MagNet and other segmentation refinement methods on the Cityscapes dataset. The backbone HRNetV2-W18+OCR <cit.> was used as the segmentation module for all refinement methods.

to @X[2.0,l]X[1,c]X[1,c]X[1,c]X[1,c]X[1,c]@
fence     pole     traffic sign     person     motor-cycle

Patching     33.32     42.87     59.39     61.23     22.12

Downsample     45.38     35.58     54.97     64.64     36.16

DenseCRF <cit.>     45.30     32.27     54.69     64.32     36.20

DGF <cit.>     45.42     35.63     55.27     64.85     36.20

PointRend <cit.>     45.01     42.71     60.18     67.10     39.17

SegFix <cit.>     46.17     38.77     59.23     67.86     37.12

MagNet-Fast     48.57     46.38     64.84     69.65     41.98

MagNet     50.59     49.39     64.15     72.16     45.19

Improvement     4.42     6.52     4.66     4.30     6.02

IoU for some specific categories on Cityscapes. The best and previous best method is highlighted in red and blue color respectively, and the difference between them is shown in the last row.

tab:cityscapes_result and tab:cityscapes_details compare the performance of MagNet with several state-of-the-art semantic segmentation refinement methods. All methods were trained and applied to the output of the pretrained HRNet+OCR <cit.> model, which is among the leading methods on this dataset. Although there are various HRNet models <cit.>, we used HRNetV2-W18 given its manageable complexity that was required for our experiments, especially in terms of memory constraint.

MagNet-Fast is an efficient version of MagNet, in which only the most uncertain patches are refined at each scale. In this experiment, we selected the number of patches so that MagNet-Fast had a similar processing speed as SegFix <cit.> and DGF <cit.>. This model ran on only three scales $256{\rightarrow}512{\rightarrow}2048$ and the three highest uncertain patches for each scale. In total, MagNet-Fast needed to run inference on $1 + 3 + 3 = 7$ patches, comparing to $1 + 4 + 16 + 64 = 85$ patches of MagNet, and 64 patches of the patch processing approach.

In this experiment, except for DenseCRF <cit.>, we finetuned other frameworks to achieve the best result with the segmentation backbone. For SegFix <cit.>, with the offset prediction published by the authors, the best result was achieved with the offset width of 10. Both DGF <cit.> and PointRend <cit.> were trained with the output of the segmentation backbone. Besides, DGF ran on patches to be fairly compared with our method in speed and memory.

As can be observed, DenseCRF <cit.> is the only method that cannot improve the coarse segmentation. PointRend <cit.> is the fastest method, but the improvement is small. The inference times of MagNet-Fast, SegFix <cit.>, and DGF <cit.> are similar, but MagNet-Fast outperforms the others significantly. MagNet was slower, but it had the highest mIoU.

The cumulative IoU distributions of these methods in our experiment are shown on fig:cityscapes_result:dist. There is a big gap between MagNet and the other methods, especially when looking at the zoom-in window.

fig:cityscapes_result:vis shows the results of several methods on two Cityscapes images. Both MagNet and MagNet-Fast yield the best refinement. SegFix <cit.> cannot recover small objects, such as sign poles, that are wrongly merged with a bigger region, while PointRend <cit.> performs poorly due to the lack of global context.

§.§.§ Ablation study: point selection

to @X[1,c]X[1,c]X[1,c]X[1,c]@
$Y^u$     $(1 - R^u)$     $\mF$     mIoU (%)

Performance of MagNet on Cityscapes with different ranking scores. The initial segmentation has mIoU of $63.23\%$. With $k{=}2^{16}$, the framework achieves the best performance when using both $Y^u$ and $(1 - R^u)$ with smoothing operation.

Correlation between the number of selective points of each scale and mean IoU on the Cityscapes dataset. When the quantity of points increases, the performance of MagNet continuously grows while the mIoU of PointRend decreases.

tab:cityscapes_ranking_score shows our ablation study on the importance of using the prediction uncertainty maps $Y^u$, $R^u$, and the median filtering function $\mF$ in selecting points for refinement. The best performance was achieved when both uncertainty maps were used. Also, smoothing with median filtering improved the result in every case.

We also studied how the number of refinement points correlates with accuracy. The results of MagNet and PointRend <cit.> are shown in fig:point_comparison. As can be seen, the performance of MagNet improved when the number of points increased. The performance stopped increasing after $2^{16}$ points, and it dropped to 66.86% if all points were selected for refinement. Meanwhile, the performance of PointRend decreased significantly when the number of selected points increased beyond 2048; it even dipped below the initial value where no refinement was applied.

§.§.§ Ablation study: segmentation backbones

to @X[3,l]X[1,c]@
Model     mIoU(%)

2lBackbone: DeepLabV3+ <cit.>

Patch processing     59.64

Downsampling     52.01

MagNet     61.99

2lBackbone: HRNetV2-W48 + OCR <cit.>

Patch processing     54.30

Downsampling     63.92

MagNet     68.90

Results of using MagNet with two backbone networks. MagNet can be used with different segmentation backbones, and improve their segmentation results. See tab:cityscapes_result for the results for using MagNet with HRNetV2-W18 backbone.

We also tested the MagNet framework with two different segmentation backbone networks, and the results are shown in tab:cityscapes_backbones. In both cases, MagNet improved the segmentation results of the original networks significantly, from 2% to 5%. In this experiment, we used four scale levels: $256{\rightarrow}512{\rightarrow}1024{\rightarrow}2048$ and the number of refinement locations for each patch was $k=2^{16}$.

§.§ DeepGlobe

to @X[4,l]X[3,l]X[2.5,l]X[1.5,r]@
Model     mIoU(%)     Mem.(MB)     Time(s)

U-net<cit.>     50.11     1813

FCN-8s<cit.>     52.86     10569

SegNet<cit.>     60.93     2645

DeepLabv3+<cit.>     63.50     1541

FPN<cit.>     67.86     1247     0.01

3lPatch processing

U-net<cit.>     46.53     1813

FCN-8s<cit.>     62.43     10569

SegNet<cit.>     68.40     2645

DeepLabv3+<cit.>     69.69     1541

FPN<cit.>     70.98     1247     0.31

DenseCRF<cit.>     70.36 0.62     1247     39.68

DGF<cit.>     70.38 0.6     1435 188     0.25

GLNet<cit.>     71.60 0.62      1865 618     0.37

PointRend<cit.>     71.78 0.8     1593 346     0.16

MagNet-Fast     71.85 0.87     1559 312     0.29

MagNet     72.96 1.98      1559 312     1.19

Segmentation results on the DeepGlobe dataset. We used the same segmentation backbone (FPN) for all refinement methods in the last part.

DeepGlobe is a dataset of high-resolution satellite images. The dataset contains 803 images, annotated with seven landscape classes, including the unknown class. Following the evaluation protocol of <cit.>, the unknown class is ignored in mIoU calculation, so there are only six classes to consider. The size of the images is $2448{\times}2448$ pixels. We used the same train/validation/test split as reported in <cit.>, with 455, 207, and 142 images for training, validation, and testing, respectively.

The Feature Pyramid Network (FPN) <cit.> with Resnet-50 backbone was used as the segmentation network as in the previous work GLNet <cit.>. We also used the same input size $508{\times}508$ as GLNet. We used three refinement stages with three scales $612{\rightarrow}1224{\rightarrow}2448$ and selected $2^{16}$ points for refinement at each scale. The results are shown in tab:deepglobe_result.
For MagNet-Fast, at each of the three scale levels, we selected the top three patches with the highest level of prediction uncertainty for refinement. PointRend <cit.> was also trained with the same segmentation backbone and it achieved higher accuracy than GLNet. Both MagNet and MagNet-Fast outperformed other methods.

§.§ Gleason

to @X[3.5,l]X[3,l]X[2.5,l]X[1.5,r]@
Model     mIoU(%)     Mem.(MB)     Time(s)

Experts     65.48     -     -

Patching     46.56     1903     2.42

Downsampling     68.90     1903     0.02

DenseCRF<cit.>     69.46 0.56     1903     141.79

DGF<cit.>     68.91 0.01     2223 320     0.29

PointRend<cit.>     68.97 0.07     2655 752     0.21

MagNet-Fast     69.75 0.85     2621 718     0.33

MagNet     70.60 1.7     2621 718     2.74

Performance of MagNet and other frameworks on Gleason dataset with PSPNet <cit.> as the backbone.

Gleason <cit.> is a medical dataset with histopathological images of prostate cancer. The task is to segment and grade lesions on ultra-high-resolution images. There are four classes in the dataset that need to be segmented: benign, Grade 3, Grade 4, and Grade 5. The dataset contains 244 images with a size of $5000{\times}5000$ pixels with segmentation labels provided by six clinical experts. The combined final label is based on majority voting. We randomly split the dataset into 195 training and 49 testing images.

PSPNet <cit.>, the highest-ranked method on the leaderboard for Gleason, was used as the segmentation network with the backbone Resnet-101. We used the input size of $512{\times}512$, and four refinement stages with four scales: $625{\rightarrow}1250{\rightarrow}2500{\rightarrow}5000$. MagNet-Fast was also run on four scales, but only on three patches with the highest level of prediction uncertainty at each scale. There are $2^{16}$ refinement points for each scale. The results of MagNet and MagNet-Fast, together with the result of the winning solution PSPNet and the mIoU agreement between medical experts, are shown in tab:gleason_result. MagNet was run with the PSPNet segmentation backbone, and it improved the performance of PSPNet by 1.7%.

§ CONCLUSIONS

We have proposed MagNet, a multi-scale segmentation framework for high-resolution images. MagNet can generate high-resolution segmentation output without exploding the GPU memory usage by dividing input images into patches. To avoid the problem of being too global or local, patches of multiple scales are considered, from the coarsest to the finest levels. MagNet has multiple segmentation stages, where the output of one stage will be used as the input for the next stage, and the segmentation output will be progressively refined. We have demonstrated the benefits of MagNet on three challenging high-resolution image datasets, where MagNet outperforms the previous state-of-the-art methods by a margin of 1% to 2% in terms of mean Intersection over Union (mIoU).