# Addendum to “Power-law decay of the fraction of the mixed eigenstates in
kicked top model with mixed-type classical phase space”

Hua Yan<EMAIL_ADDRESS>CAMTP - Center for Applied Mathematics and
Theoretical Physics, University of Maribor, Mladinska 3, SI-2000 Maribor,
Slovenia  Qian Wang  CAMTP - Center for Applied Mathematics and Theoretical
Physics, University of Maribor, Mladinska 3, SI-2000 Maribor, Slovenia
Department of Physics, Zhejiang Normal University, Jinhua 321004, China  Marko
Robnik<EMAIL_ADDRESS>CAMTP - Center for Applied Mathematics and
Theoretical Physics, University of Maribor, Mladinska 3, SI-2000 Maribor,
Slovenia

###### Abstract

By using the Krylov subspace technique to generate the spin coherent states in
kicked top model, a prototype model for studying quantum chaos, the accessible
system size for studying the Husimi functions of eigenstates can be much
larger than that reported in the literature and our previous study Phys. Rev.
E 108, 054217 (2023). In the fully chaotic kicked top, we find that the mean
Wehrl entropy localization measure approaches the prediction given by the
Circular Unitary Ensemble. In the mixed-type case, we identify mixed
eigenstates by the overlap of the Husimi function with regular and chaotic
regions in classical compact phase space. Numerically, we show that the
fraction of mixed eigenstates scales as $j^{-\zeta}$, a power-law decay as the
system size $j$ increases, across nearly two orders of magnitude. This
provides supporting evidence for the principle of uniform semiclassical
condensation of Husimi functions and the Berry-Robnik picture in the
semiclassical limit.

##  I Introduction

At the very beginning of quantum chaos it was Percival [1] who proposed to
distinguish between the regular quantum eigenstates and the chaotic ones,
associated with classical regular and chaotic regions in the phase space,
respectively. They should be distinguished by the nature of their energy
spectra, namely by their sensitivity with respect to small perturbations. This
line of thought led to the semiclassical eigenfunction hypothesis by Berry
[2], which is crucial for understanding quantum eigenstates. It proposes that
in the semiclassical limit, quantum eigenstates concentrate on compact phase
space regions explored by typical trajectories in the long time limit. In
integrable systems, these are the invariant tori, while in ergodic dynamics,
eigenstates equidistribute on the energy surface. In mixed-type systems with
both regular and chaotic motions in phase space [3], eigenstates concentrate
either on regular or chaotic regions. This hypothesis finally develops into
the so-called principle of uniform semiclassical condensation (PUSC) of Wigner
functions (or Husimi functions) [4]. This principle, in turn, directly leads
to the eigenstate thermalization hypothesis conjecture by Srednicki [5], which
is now recognized as an explanation for the emergence of thermodynamic
equilibrium in isolated quantum many-body systems [6, 7, 8]. Moreover, as a
consequence of this principle, the spectral statistics for regular and chaotic
states in the semiclassical limit are separately described by Poissonian
statistics [9] and random matrix theory [10, 11, 12], while the entire
spectrum is effectively captured by the Berry-Robnik picture [13].

Regarding the behavior of quantum eigenstates in mixed-type systems away from
the semiclassical limit, one expects various tunneling processes between
different phase-space structures. Among these are the chaos-assisted tunneling
between different regular regions [14, 15, 16], the flooding of chaotic states
into regular islands [17], and the so-called resonance-assisted tunneling
between the regular island across the chaotic sea [18, 19, 20, 21]. The
general behavior of quantum states approaching the semiclassical limit is not
well studied [22]. Through the overlap of the Husimi function with both
regular and chaotic regions in classical compact phase space, mixed
eigenstates can be identified, distinct from regular and fully chaotic states.
Very recently, in quantum billiards [23], it was found that the fraction of
mixed states shows a power-law decay as the semiclassical limit is approached,
across two orders of magnitude. In quantum kicked top (QKT), a prototype
Floquet model for studying quantum chaos [24, 25], our previous study [26]
also showed a power-law decay for mixed-type QKT, but only within a relatively
small range of system sizes. This limitation arises from our method used to
generate the spin coherent state (SCS), where a direct calculation of large
factorials are unavoidable.

In this addendum paper, we use the Krylov subspace techniques [27, 28, 29, 30]
to generate the SCS and avoid the factorials. This allows for a significantly
larger accessible system size for studying Husimi functions, compared to that
reported in literature and our previous study. In the fully chaotic QKT, we
find that the mean Wehrl entropy localization measure (ELM) approaches the
prediction given by the Circular Unitary Ensemble (CUE). In the mixed-type
case, we show a power-law decrease of the fraction of mixed eigenstates as
system size increases over nearly two orders of magnitude, corroborating the
PUSC and the Berry-Robnik picture in the semiclassical limit.

The paper is organized as follows. In Sec. II, we introduce QKT and analyze
the transition to chaos for both its classical correspondence and quantum
case. In Sec. III, we give the definition of the SCS and apply the Krylov
subspace technique to generate them. We then establish a consistent transition
to chaos using the mean Wehrl ELM, compared to other indicators of chaos. In
the fully chaotic QKT, we study the mean Wehrl ELM versus the system size. In
Sec. IV, we analyze the joint distribution of the phase-space overlap index
and Wehrl ELM, and show as the system size increases, how the distribution of
the former evolves, as well as the fraction of mixed states. We draw our
conclusion in Sec. V.

##  II Quantum kicked top

The QKT is described by the Hamiltonian [24]

$\displaystyle H=\alpha
J_{x}+\frac{\gamma}{2j}J_{z}^{2}\sum_{n=-\infty}^{\infty}\delta(t-n),$ (1)

where the dynamical variables of the top are three components of the angular
momentum operators of the spin-$j$ system, can also be expressed in terms of
$2j$ collective spin-$\frac{1}{2}$ Pauli operators, for example
$J_{z}=\sum_{k=1}^{2j}\sigma_{z}^{(k)}/2$. The dimension of the Hilbert space
is $N=2j+1$, and the squared angular momentum is conserved, $J^{2}=j(j+1)$
with $j$ integer or half-integer. The first term in Eq. (1) describes a
precessional rotation about the $x$-axis with angular frequency $\alpha$, the
second term denotes a torsional rotation around the $z$-axis with strength
$\gamma$ (with the period set to unity).

The dynamical evolution of the QKT is governed by the Floquet operator

$\displaystyle F=\exp(-i\frac{\gamma}{2j}J_{z}^{2})\exp(-i\alpha J_{x}),$ (2)

is invariant under rotation by the angle $\pi$ along the $x$-axis
$R_{x}=e^{i\pi J_{x}}$, which suggests the even or odd parity of the quasi-
energy. Using this symmetry, the $(2j+1)$-dimensional unitary matrix $F$ in
the representation of the Dicke states {$|j,m\rangle,|m|\leq j$}, i.e.
eigenstates of $J_{z}$, can further be reduced to two unitary matrices

$\displaystyle F_{m,m^{\prime}}^{\pm}$ $\displaystyle=\langle
j,m,\pm|F|j,m^{\prime},\pm\rangle,$ (3)

where $|j,m,\pm\rangle=(|j,m\rangle\pm|j,-m\rangle)/\sqrt{2}$, $1\leq m\leq j$
for the odd parity, and $0\leq m\leq j$ for the even parity where
$|j,0,+\rangle$ is not normalized. It should be noted that at $\alpha=\pi/2$,
the QKT has a further symmetry [25]. In this addendum paper, we set
$\alpha=11\pi/19$, the same as the previous work.

The Heisenberg equation of motion is given by the map
$J_{i,n+1}=F^{\dagger}J_{i,n}F$, where $J_{i,n}\equiv J_{i}(n)$ denotes the
time-evolved operators at $t=n$. Using Baker-Campbell-Hausdorff (BCH) formula
for the expansion of the map, one has

$\displaystyle J_{x,n+1}$
$\displaystyle=\frac{1}{2}[\tilde{J}_{x,n}+i\tilde{J}_{y,n}]\exp\left[i\frac{\gamma}{2j}(2\tilde{J}_{z,n}+1)\right]+h.c.,$
$\displaystyle J_{y,n+1}$
$\displaystyle=\frac{1}{2i}[\tilde{J}_{x,n}+i\tilde{J}_{y,n}]\exp\left[i\frac{\gamma}{2j}(2\tilde{J}_{z,n}+1)\right]+h.c,$
$\displaystyle J_{z,n+1}$ $\displaystyle=\tilde{J}_{z,n},$ (4)

where $\tilde{J}_{x,n}=J_{x,n}$, $\tilde{J}_{y,n}=J_{y,n}\cos\alpha-
J_{z,n}\sin\alpha$, and $\tilde{J}_{z,n}=J_{y,n}\sin\alpha+J_{z,n}\cos\alpha$.

###  II.1 The transition to Hamiltonian chaos

The classical map emerges in the classical limit as $j\to\infty$ where
noncummutivity is negligible and the rescaled Heisenberg operators
$\mathbf{x}_{n}\equiv\mathbf{x}(n)=\mathbf{J}_{n}/j$ tend to a unit vector
with commuting components. From Eq. (II), the classical equations of motion
follow

$\displaystyle x_{n+1}$
$\displaystyle=\tilde{x}_{n}\cos\gamma\tilde{z}_{n}-\tilde{y}_{n}\sin\gamma\tilde{z}_{n},$
$\displaystyle y_{n+1}$
$\displaystyle=\tilde{x}_{n}\sin\gamma\tilde{z}_{n}+\tilde{y}_{n}\cos\gamma\tilde{z}_{n},$
$\displaystyle z_{n+1}$ $\displaystyle=\tilde{z}_{n},$ (5)

where $\tilde{x}_{n}=x_{n}$, $\tilde{y}_{n}=y_{n}\cos\alpha-z_{n}\sin\alpha$,
and $\tilde{z}_{n}=y_{n}\sin\alpha+z_{n}\cos\alpha$. This classical map
$\mathbf{x}_{n+1}=f(\mathbf{x}_{n})$ is restricted on the unit sphere
$\mathbb{S}^{2}$ and the Jacobian matrix $\partial
f/\partial\mathbf{x}_{n}=\Omega_{n}\cdot R_{\alpha}$, where

$\displaystyle\Omega_{n}$
$\displaystyle=\begin{pmatrix}\cos\gamma\tilde{z}_{n},&-\sin\gamma\tilde{z}_{n},&-\gamma\tilde{x}_{n}\sin\gamma\tilde{z}_{n}-\gamma\tilde{y}_{n}\cos\gamma\tilde{z}_{n}\\\
\sin\gamma\tilde{z}_{n},&\cos\gamma\tilde{z}_{n},&\gamma\tilde{x}_{n}\cos\gamma\tilde{z}_{n}-\gamma\tilde{y}_{n}\sin\gamma\tilde{z}_{n}\\\
0&0&1\end{pmatrix},$ $\displaystyle R_{\alpha}$
$\displaystyle=\begin{pmatrix}1&0&0\\\ 0&\cos\alpha&-\sin\alpha\\\
0&\sin\alpha&\cos\alpha\end{pmatrix}.$ (6)

Figure 1: (a) A stroboscopic map of the classical dynamics generated from 20
random initial conditions of $5\times 10^{4}$ iterations (top panels). (b) The
corresponding logarithmic values of SALI (bottom panels) on the parametrized
phase space $(\theta,\phi)$ discretized by $200\times 400$ grids of square
cells, of same area, after 300 kicks. The initial conditions colored dark blue
correspond to chaotic orbits, the yellowish indicates ordered motion, and the
intermediate suggests sticky orbits. From left to right, the kicking strength
$\gamma=2,4,6$ and $\alpha=11\pi/19$.

The stroboscopic dynamics of the classical kicked top (CKT) is symplectic
since $\det(\partial f/\partial\mathbf{x}_{n})=1$, and $\mathbb{S}^{2}$ as the
phase space of the classical top can be parametrized as:
$x=\sin\theta\cos\phi$, $y=\sin\theta\sin\phi$, $z=\cos\theta$, with $z$ and
$\phi$ being the canonical pair. This map undergoes a transition from
integrability to chaos with increasing kicking strength $\gamma$, as
illustrated in Fig. 1(a) the Poincaré section of the classical map, for
different values of $\gamma$. In the classical phase space, regular orbits
predominate at small values of $\gamma$. As the kicking strength increases,
mixed dynamics emerge, characterized by the coexistence of regular orbits and
chaotic motion. Eventually, for larger $\gamma$, the system transitions into
full chaos, as illustrated in the rightmost panel of Fig. 1(a).

To qualitatively demonstrate the transition to chaos in the classical map, in
this addendum, we use the smaller alignment index (SALI) [31, 32, 33] which
rely on the evolution of deviation vectors from a given orbit, to detect
regular and chaotic motion of the symplectic map. For chaotic orbits, it was
proved that SALI$(t)\propto e^{(-L_{1}-L_{2})t}$ [34], with $L_{1},L_{2}$
being the two largest Lyapunov exponents (LEs). While in three-dimensional
symplectic map $L_{2}=0$, therefore for CKT, the SALI is directly related to
the largest LE as SALI$(t)\propto e^{-L_{1}t}$. Fig. 1(b) illustrates the
logarithmic SALI values after 300 kicks for grid points, at each point plotted
using an assigned color according to its value. Upon comparing with the
corresponding Poincaré sections in Fig. 1(a), in addition to the resemblance,
we also acquire a much more detailed picture of the regions where chaotic or
regular motion occurs, and at the borders between these regions we find
intermediate colors which correspond to sticky orbits.

To further quantify the degree of chaos, we define the chaotic fraction
$\mu_{c}$ as the relative area of the chaotic part of the phase-space as

$\displaystyle\mu_{c}=\frac{1}{4\pi}\int\chi_{c}(\theta,\phi)\sin\theta
d\theta d\phi,$ (7)

where $\chi_{c}(\theta,\phi)$ denotes the characteristic function of the
chaotic component, which takes the value of 1 on chaotic region and zero
otherwise. The chaotic fraction $\mu_{c}$ can be used as an indicator of
chaos. It measures the transition from integrable dynamics with $\mu_{c}$ = 0
to the fully chaotic (ergodic) $\mu_{c}=1$. Our criterion for the
classification of initial conditions belonging to a chaotic region is that
SALI$(t)\leq 10^{-8}$ at $t=300$, implying that the deviation vectors have
been aligned. Calculated from $200\times 400$ initial grid points on the
parametrized phase space $(\theta,\phi)$, Fig. 2(a) illustrates the dependence
of the chaotic fraction $\mu_{c}$ on $\gamma$, for CKT. The value of $\mu_{c}$
remains near 0 for $\gamma\leq 2$ and increases thereafter until it saturates
to 1 for $\gamma\gtrsim 4.4$. It should be noted that the transition point at
$\gamma\simeq 2$ obtained by the chaotic fraction $\mu_{c}$ from the SALI
method is consistent with the transition point from the Kolmogorov-Sinai
entropy given in Ref. [26].

Figure 2: (a) Chaotic fraction $\mu_{c}$ of the classical phase space as a
function of $\gamma$, calculated from $200\times 400$ initial grid points on
the parametrized phase space $(\theta,\phi)$. (b) The normalized mean spacing
ratio $r_{c}$ averaged over the odd and even parity, versus $\gamma$ for QKT
with $j=2^{12}$. In both figures, the lower and upper (red) dashed lines
denote the value $0$ and $1$, and the parameter $\alpha=11\pi/19$.

###  II.2 The transition to quantum chaos

The transition to Hamiltonian chaos, as indicated by the variation of
$\mu_{c}$ versus the kicking strength in Fig. 2(a), would also be manifested
in the quantum realm due to the quantum-classical correspondence. There are
several methods to diagnose quantum chaos, statistically or dynamically. One
is the short-range statistical properties of spectra, particularly the
nearest-neighbor-level spacing distribution, denoted as $P(s)$, where $s$
denotes the spacing between two unfolded consecutive energy levels. In the
quantum integrable case, this distribution follows the Poisson distribution
$P(s)=e^{-s}$, whereas the quantum chaotic case with the time-reversal
symmetry is the well-known Wigner-Dyson distribution $P(s)=(\pi/2)s\exp(-\pi
s^{2}/4)$ from the random-matrix statistics.

To avoid the numerical spectral unfolding, we consider instead the spacing
ratios, defined as [35, 36]

$\displaystyle r_{n}=\min(\frac{s_{n}}{s_{n-1}},\frac{s_{n-1}}{s_{n}}),$ (8)

with $s_{n}=\nu_{n}-\nu_{n-1}$ the consecutive level spacing, from
$\\{\nu_{n}\\}$ an ordered set of quasienergy levels, that
$F|\nu_{n}\rangle=e^{i\nu_{n}}|\nu_{n}\rangle$. The approximate formulas for
the distribution of $P(r)$ can be derived from random matrix statistics. For
the Poisson level spacing, the resulting mean spacing ratio $\langle r\rangle$
is $\langle r\rangle_{p}=2\ln 2-1\approx 0.386$, whereas for Circular
Orthogonal Ensemble (COE) statistics, $\langle r\rangle_{COE}\approx 0.527$
[37]. We can then define a quantity $r_{c}$ the normalized mean spacing ratio
as an indicator of quantum chaos, namely

$\displaystyle r_{c}=\frac{\langle r\rangle-\langle r\rangle_{P}}{\langle
r\rangle_{COE}-\langle r\rangle_{P}}.$ (9)

Clearly, $r_{c}=0$ corresponds to the integrable case, while $r_{c}=1$
signifies full chaos. As discussed in Ref. [25], in QKT, there exist two
antiunitary time reversal operators

$\displaystyle T_{1}=e^{i\alpha J_{x}}e^{i\pi J_{z}}K,\quad T_{2}=e^{-i\alpha
J_{x}}e^{i\pi J_{y}}K,$ (10)

where $K$ denotes the conjugation in the eigenbasis of $J_{z}$. These
operators satisfy $T_{1}^{2}=T_{2}^{2}=1$ and time reversal property
$T_{1}FT_{1}=T_{2}FT_{2}=F^{\dagger}$. One can therefore expect that in the
fully chaotic regime, the statistics of quasienergy of $F$ must be given by
COE. As depicted in Fig. 2(b), the variation of $r_{c}$, averaged from both
parities, with kicking strength $\gamma$ illustrates that it remains close to
0 for $\gamma\leq 2$, gradually increasing thereafter until it reaches a
plateau near 1 for $\gamma\geq 4.4$. This behavior precisely aligns with the
chaotic fraction $\mu_{c}$ versus $\gamma$ depicted in Fig. 2(a), as well as
with the mean Wehrl entropy localization length of eigenstates, which we will
demonstrate in the following.

##  III Husimi function and Localization of eigenstates

By employing spectral statistics, we have demonstrated in Sec. II.2 the
transition to quantum chaos, which aligns with the classical route to chaos.
In this section, we introduce the SCS and employ the Krylov subspace technique
to numerically calculate the SCS. Notably, by avoiding the problem of large
factorials, the accessible value of $j$ can be several orders of magnitude
higher than those reported in the literature, including our previous work
[26]. Following this, by analyzing the Husimi function of eigenstates based on
SCS, we can explore the localization of these eigenstates, especially the
Wehrl entropy [38] and the corresponding entropy localization measure (ELM),
to a deeper semiclassical regime.

###  III.1 Spin coherent state and Husimi functions

The spin coherent state $|\theta,\phi\rangle$ from SU(2) [39, 40, 41], can be
obtained as a rotation on the Dicke manifold, expanded in the Dicke basis
$|j,m\rangle$ as

$\displaystyle|\theta$
$\displaystyle,\phi\rangle=R(\theta,\phi)|j,j\rangle=\exp(\mu
J_{-}-\mu^{*}J_{+})|j,j\rangle$ (11a)
$\displaystyle=(1+|\xi|^{2})^{-j}\sum_{m=-j}^{j}\binom{2j}{j+m}^{1/2}\xi^{j-m}|j,m\rangle,$
(11b)

where $\mu=\frac{\theta}{2}\exp(i\phi)$,
$\xi=\tan(\frac{\theta}{2})\exp(i\phi)$ with $\theta\in[0,\pi]$ and
$\phi\in[-\pi,\pi)$, $J_{\pm}=J_{x}\pm iJ_{y}$. The SU(2) generalized coherent
states $|\theta,\phi\rangle$ are overcomplete that

$\displaystyle\frac{2j+1}{4\pi}\int\sin\theta d\theta d\phi\
|\theta,\phi\rangle\langle\theta,\phi|=\mathbb{I}.$ (12)

The Husimi function of the eigenstate $|\nu_{n}\rangle$ of the Floquet
operator $F$ is then given by

$\displaystyle Q_{n}(\theta,\phi)=|\langle\theta,\phi|\nu_{n}\rangle|^{2},$
(13)

with the normalization condition

$\displaystyle\frac{2j+1}{4\pi}\int\sin\theta d\theta d\phi\
Q_{n}(\theta,\phi)=1.$ (14)

Numerical studies of Husimi functions in the past usually employed Eq. (11b)
for the SCS representation of quantum states. However, the presence of large
factorials in this equation imposes limitations on the attainable values of
$j$, typically only a few hundred. A further examination of Eq. (11a) shows
that the SCS can be obtained directly from a unitary transformation of the
polar state $|\theta=0,\phi\rangle=|j,j\rangle$. This unitary transformation
is an exponential of a sparse matrix, because in the Dicke basis

$\displaystyle J_{\pm}|j,m\rangle=\sqrt{(j\mp m)(j\pm m+1)}|j,m\pm 1\rangle.$
(15)

Exponential of sparse matrix on vectors can be effectively computed using
Krylov subspace technique (see Ref. [28]), this alternative approach bypasses
the problem of large factorials and enables the attainment of significantly
larger value of $j$ (differing by several orders of magnitude, depending on
the sparsity), compared with the previous studies. One can take an initial
look at the Husimi functions of eigenstates provided in Sec. IV.1 and more
illustrations in the Appendix, where $j=2^{12}$.

###  III.2 Localization of eigenstates

One common way for defining the localization of eigenstates is via Shannon
entropy. The eigenstates from odd parity $|\nu_{n}\rangle_{odd}$ that are
expressed in the basis of unperturbed Hamiltonian, as

$\displaystyle H_{n}=-\sum_{m=1}^{j}|\nu_{nm}|^{2}\ln|\nu_{nm}|^{2},\
\nu_{nk}=\langle x_{m}|\nu_{n}\rangle_{odd},$ (16)

where $|x_{m}\rangle$ is the eigenbasis of the unperturbed Floquet part
$\exp({-i\alpha J_{x}})$, i.e., the eigenbasis of $J_{x}$ in odd parity. For
the even parity, it can be defined in a same way, and we have verified that
the statistics of $H_{n}$ is nearly identical for both parity. This is also
referred to as the eigenvector statistics. The Shannon ELM is then defined as
$L_{n}=e^{H_{n}}/N_{odd}$ with $N_{odd}=j$, and the mean localization measure
is denoted as $\langle L\rangle$. Apart from this basis-dependent method,
another approach for measuring the localization of eigenstates is basis-
independent, relying on the Wehrl entropy of eigenstates. It is defined as

$\displaystyle S_{n}=-\frac{2j+1}{4\pi}\int\sin\theta d\theta d\phi\
Q_{n}(\theta,\phi)\ln Q_{n}(\theta,\phi).$ (17)

The Wehrl ELM is $\mathcal{L}_{n}=e^{S_{n}}/N$, with $N=2j+1$ and
$\langle{\mathcal{L}}\rangle$ denotes the averaging. It’s worth noting that
Wehrl conjectured that the minimum of Wehrl entropy is attained for coherent
state [42]. Lieb originally proved the conjecture [43], later extended it to
SU(2) and SU(N) coherent states [44, 45, 46]. However, the problem of the
uniqueness of minimizers in some cases is still open [47].

One may anticipate that a random eigenstate will tend to the uniform
distribution that $|\nu_{nk}|^{2}=1/N_{odd}$ and $Q_{n}(\theta,\phi)=1/N$,
resulting in the maximum value of both entropy $H_{n}=\ln N_{odd}$ and
$S_{n}=\ln N$, as well as ELMs $L_{n}=\mathcal{L}_{n}=1$. For a random pure
state $|\mathcal{R}\rangle$ from COE/GOE [48, 49], it was proved that mean
Shannon entropy $\langle H\rangle_{\mathcal{R}}=\Psi(N/2+1)-\Psi(3/2)\sim\ln
N+\gamma_{e}+\ln 2-2\approx\ln N-0.7296$, with $\Psi(x)$ the digamma function
and $\gamma_{e}$ denotes the Euler constant, where the approximation is for
the asymptotic limit $N\to\infty$. Additionally, the mean Shannon/Wehrl
entropy of a random pure state from CUE is $\langle
S\rangle_{\mathcal{R}}=\Psi(N+1)-\Psi(2)\sim\ln N+\gamma_{e}-1\approx\ln
N-0.4228$ [50]. The random pure state statistics then gives $\langle
L\rangle_{\mathcal{R}}\approx 0.482$, and
$\langle\mathcal{L}\rangle_{\mathcal{R}}\approx 0.655$ for large $N$.

Figure 3: Mean Shannon ELM $\langle L\rangle$ from the odd parity and mean
Wehrl ELM $\langle\mathcal{L}\rangle$, versus $\gamma$ for QKT with $j=2^{12}$
and $\alpha=11\pi/19$. The dashed horizontal lines indicate the mean Shannon
ELM of COE and Wehrl/Shannon ELM of CUE.

In Fig. 3, we study the mean Shannon ELM $\langle L\rangle$ in the basis of
$J_{x}$, for the odd parity (we found for the even parity, it is the same) and
mean Wehrl ELM $\langle\mathcal{L}\rangle$. We computed the Wehrl ELM using a
grid of $N_{p}=200\times 400$ points in the $(\theta,\phi)$ space, with each
point holding a coherent state, where all the Husimi functions are normalized
according to the discretization of Eq. (14), as

$\displaystyle\frac{(2j+1)\pi}{2N_{p}}\sum_{ij}Q_{n}(\theta_{i},\phi_{j})\sin\theta_{i}=1.$
(18)

It shows that the Shannon ELM saturates to the COE value $\langle
L\rangle_{R}$ for $\gamma\gtrsim 4.4$, and before that, it exhibits two
different types of growth: for $\gamma>2$, it grows faster than for
$\gamma\leq 2$. While the transition versus $\gamma$ is more pronounced in the
mean Wehrl ELM. For $\gamma\leq 2$, it maintains a value equal to
$\langle\mathcal{L}\rangle_{J_{x}}\sim N^{-1/2}$, representing the mean Wehrl
ELM of the operator $J_{x}$ (see numerical result shown in Fig. 9 and proof in
Appendix A), which is close to 0 for large $N$. It then increases until it
saturates for $\gamma\gtrsim 4.4$, to a value close to
$\langle\mathcal{L}\rangle_{R}$ from CUE.

Numerical computations show that the mean Wehrl entropy of COE is slightly
smaller than that $\langle S\rangle_{\mathcal{R}}$ from CUE, and it was
conjectured in Ref. [50] that the difference vanishes in the asymptotic limit
$N\to\infty$. However, both a numerical verification at larger $N$ and a
formal mathematical proof for this conjecture is lacking. In Fig. 4, using the
Krylov subspace technique generated coherent states, we study the mean Shannon
entropy $\langle H\rangle$ and Wehrl entropy $\langle S\rangle$, of the KT at
$\gamma=8$ where the system is fully chaotic. The eigenvector statistic
reveals that $\langle H\rangle\simeq\ln N-0.7322$, is nearly equal to $\langle
H\rangle_{\mathcal{R}}$ for COE. Similarly, for Shannon ELM, $\langle
L\rangle\simeq\langle L\rangle_{\mathcal{R}}$, as shown in the inset (bottom
right). On the other hand, the mean Wehrl entropy $\langle S\rangle\simeq\ln
N-0.4573$, is close to $\langle S\rangle_{\mathcal{R}}$ for CUE, while the
inset in the top left illustrates that $\langle\mathcal{L}\rangle$ approaches
the CUE value $\langle\mathcal{L}\rangle_{\mathcal{R}}$ with increasing system
size. We have examined $\langle\mathcal{L}\rangle$ against $N$ with a further
increase in grid points to $N_{p}=250\times 500$, and it exhibits a similar
trend.

Figure 4: Mean Shannon entropy $\langle H\rangle$ and mean Wehrl entropy
$\langle S\rangle$, versus $N$ for QKT, with $\gamma=8$ and $\alpha=11\pi/19$,
where for the former $N=N_{odd}=j$ and the latter $N=2j+1$. Two dashed lines
are from best fitting. Two insets show two mean ELMs versus $N$, respectively,
where the dashed lines indicate $\langle\mathcal{L}\rangle_{\mathcal{R}}$ from
CUE (top left) and $\langle L\rangle_{\mathcal{R}}$ from eigenvector
statistics of COE (bottom right). Figure 5: (a) A joint distribution
$(\mathcal{L},M)$ for 8193 eigenstates of QKT Floquet operator, with
$\gamma=2.6$, $\alpha=11\pi/19$ and $j=2^{12}$. Husimi functions of
eigenstates isolated by boxes at different positions are plotted from panel
(a1) to (d3), the darkest blues show the area where $Q_{n}<10^{-16}$. (b)
Logarithmic values of SALI for the corresponding CKT, on the parametrized
phase space $(\theta,\phi)$ after 300 kicks.

##  IV Phase-space overlap index

The Husimi functions of eigenstates in quantum maps [51, 52], are precise
analog of the classical Poincaré section. From this correspondence, the Husimi
functions of eigenstates can be used for the identification of regular, mixed,
and chaotic eigenstates, by the criterion of overlap with the classical
Poincaré section, where one can employ the SALI method to identify whether an
initial condition belongs to the classical chaotic or regular regions. As it
was introduced and implemented in previous works [53, 54, 55], for eigenstate
$|\nu_{n}\rangle$, the overlap can be quantified by the overlap index $M_{n}$,
defined as

$\displaystyle M_{n}=\frac{2j+1}{4\pi}\int\sin\theta d\theta d\phi\
f(\chi_{c})Q_{n}(\theta,\phi),$ (19)

where $f(\chi_{c})=1$ if $\chi_{c}=1$, and $f(\chi_{c})=-1$ if $\chi_{c}=0$.
Here, $\chi_{c}(\theta,\phi)$ denotes the characteristic function of the
chaotic component, as defined in Eq. (7).

According to PUSC in mixed-type systems, for Husimi functions of eigenstates,
$M$ should be either $+1$ or $-1$ in the strict semiclassical limit,
corresponding to chaotic or regular eigenstates, respectively. However, since
the semiclassical limit is not easily reached in practice, $M$ actually varies
between $-1$ and $+1$. Upon approaching the semiclassical limit, a reduction
in the fraction of mixed eigenstates with $|M|<1$ would occur, as evidenced
first in the billiard system [23], particularly exhibiting a power-law decay
with increasing energy, i.e. decreasing effective Planck constant. In the
following, we analyze quantum eigenstates of Floquet operator in mixed-type
QKT by studying the joint distribution of Wehrl ELM and overlap index. Then,
we extend the analysis of the power-law decay of the fraction of mixed
eigenstates as system size increases, as conducted in Ref. [26], up to
$j=2^{13}$. For numerical calculation of Husimi functions, we use a grid of
$N_{p}=200\times 400$ points in the $(\theta,\phi)$ space, with each point
holding a coherent state.

###  IV.1 The overlap index and Wehrl ELM

From the definition of Wehrl ELM following Eq. (17), regular eigenstates
characterized by an overlap index of $M=-1$ and condensing on invariant tori
are localized, yielding lower values of $\mathcal{L}$. In contrast, chaotic
eigenstates with $M=+1$ can be either strongly localized or extended. Previous
studies have numerically shown that the probability distribution
$P(\mathcal{L})$ of chaotic eigenstates indeed follows a beta distribution
[56, 57, 58]. As the system approaches the semiclassical limit, it tends
towards a delta distribution
$\delta(\mathcal{L}-\langle\mathcal{L}\rangle_{\mathcal{R}})$. It is important
to note that in mixed-type systems, to verify this distribution, the Wehrl ELM
of the chaotic eigenstates should be further renormalized by the relative area
$\mu_{c}$ of the chaotic component in the phase space. In addition to strictly
regular and chaotic states, there would exist mixed eigenstates with $|M|<1$
in mixed-type systems, exhibiting varying degrees of localization. This is
illustrated in Fig. 5(a) the joint distribution $(\mathcal{L},M)$ for KT at
$\gamma=2.6$, alongside selected representative Husimi functions of
eigenstates from different regions of the parameter space, as shown in Fig.
5(a1)-(d3). Moreover, Fig. 5(b) presents the SALI plot from the corresponding
classical dynamics at $\gamma=2.6$ in phase space, for comparison. At the
border of the regular islands and chaotic sea, visible stickiness are denoted
by intermediate colors.

Figure 6: Histograms of $P(M)$ for increasing system size, from left to right.
In panels (a1)-(a4), the kicking strength is set to $\gamma=2.3$. In panels
(b1)-(b4), it is $\gamma=2.6$, and in panels (c1)-(c4), it is $\gamma=3$.

Fig. 5 shows that states predominantly cluster at both ends of the overlap
index, one end representing regular states and the other chaotic states, while
mixed states sparsely distribute in between. This clustering is more evident
for QKT at $\gamma=3$, as shown in Fig. 11 (see Appendix B the gallery of
states), where larger values of $\gamma$ correspond to increased $\mu_{c}$ the
fraction of chaotic components, and a simpler structure of classical phase
space. Examining Husimi functions of the chaotic eigenstates plotted in Fig.
5(a1)-(a3), alongside the corresponding SALI plot, a reduction in the area of
the darkest blue regions (with $Q_{n}\leq 10^{-16}$) belonging to regular
islands can be observed. The most extended state shown in Fig. 5(a3),
infiltrates into the regular islands due to quantum tunneling, and exhibits a
minor overlap with the outer tori of the regular islands. The state with the
maximum overlap index, as plotted in Fig. 5(a1), is strongly localized around
and avoids the outer tori of the regular region. Similar observation extends
to another layer of overlap index $M\simeq 0.5$ for states illustrated in Fig.
5(b1)-(b3), and to $M\simeq 0$ in Fig. 5(c1)-(c3), except that the most
localized ones now live in the vicinity of the regular islands, with smaller
contribution from the chaotic sea. Fig. 5(d1)-(d3) illustrate states from the
bottom layer of overlap index. Particularly, Fig. 5(d1) shows a state
predominantly localized at a chain of regular islands, which exhibits
tunneling between these regular islands. In Fig. 5(d2), the state also resides
in the vicinity of the regular island, with visible stickiness surrounding it,
appearing more extended. The most localized state shown in Fig. 5(d3) is
located exactly at the fixed points.

In summary, the joint distribution $(\mathcal{L},M)$ offers a comprehensive
overview of mixed-type QKT eigenstates. The study of Husimi functions shows
that, mixed states predominantly occupy chaotic regions with $M>0$, flooding
into regular areas, while states concentrate in regular islands with $M<0$,
displaying tunneling between them. Furthermore, among states with a smaller
overlap index, the Wehrl ELM shows a narrower range.

###  IV.2 The fraction of mixed eigenstates

In the semiclassical limit, it is expected that mixed states gradually
disappear in accordance with the Berry-Robnik picture and PUSC. This is
supported by the histograms depicted in Fig. 6, which illustrate the
distribution $P(M)$ of the overlap index $M$ as the system size $j$ increases,
at different kicking strengths. As the semiclassical limit is approached with
increasing $j$, we observe that: (i) Increasing values of $j$ lead to a more
pronounced alignment of eigenstates with either regular or chaotic clusters.
Specifically, at $\gamma=2.3$, the kicking strength closer to the
integrability breaking at $\gamma\simeq 2$, the chaotic cluster centers around
$M\simeq 0.8$, correlating with a much more complex phase space structure (see
Fig. 10 in Appendix B) and the possible classical partial transport barriers
[59, 60, 61, 62, 63, 64]. (ii) Fluctuations among intermediate values of $M$
decrease as the system size $j$ increases. (iii) Increasing $\gamma$ enhances
chaos and $\mu_{c}$, reflected in $P(M)$ with a higher peak at $M=+1$ and a
lower peak at $M=-1$, and reduced fluctuations for $|M|<1$.

Figure 7: Decay of the fraction of mixed eigenstates $\chi_{M}$ with respect
to $j$, at three kicking strengths $\gamma=2.3$ (squares), $\gamma=2.6$
(circles) and $\gamma=3$ (triangles). The mixed eigenstates criteria here is
$M\in[-0.9,0.5]$. The dashed lines show the power law decay of $\chi_{M}\sim
j^{-\zeta}$.

Following these observations, to quantify the decay of the fraction of mixed
states with respect to the increasing system size $j$, we define

$\displaystyle\chi_{M}=N(M\in[M_{0},M_{1}])/N,$ (20)

as the fraction of eigenstates within an interval of the overlap index
$M_{0}\leq M\leq M_{1}$, given the total number of states $N=2j+1$. In Fig. 7,
we show the case where $M\in[-0.9,0.5]$ and plot $\chi_{M}$ as a function of
$j$ for different $\gamma$. This behavior is well described (fitted) by a
power-law decay, $\chi_{M}\sim j^{-\zeta}$, with the power exponent $\zeta$
being very similar for different $\gamma$. Here, the range of $j$ we
considered is expanded by two orders of magnitude compared to the previous
study [26], which only examined $j$ up to 400, resulting in substantially
improved fitting.

The decay exponent $\zeta$ certainly depends on the interval we consider for
$M$. The interval should not be too small, lacking sufficient states for good
statistics, or too large to show little variation across different intervals.
In Fig. 8, we illustrate the variation of the decay exponent when considering
small intervals $M\in[M-\Delta M/2,M+\Delta M/2]$, with $\Delta M=0.3$, for
different $\gamma$. The inset displays the goodness of the best-fitting,
indicated by $R^{2}$, the coefficient of determination. Disregarding a few
imperfect fittings, the data suggests a decay exponent range of
$0.2<\zeta<0.45$. Given the similar power-law decay observed in the fraction
of mixed eigenstates in both billiards [23] and quantum Henon-Heilés system
[58], along with the Dicke model [65], we conjecture that this behavior is a
universal property of mixed states, as the system approaches the semiclassical
limit. A more profound open question is whether there exists an upper bound of
the decay exponent.

Figure 8: Power-law decay exponent $\zeta$ for $M\in[M-\Delta M/2,M+\Delta
M/2]$, with $\Delta M=0.3$, across different intervals. The inset shows the
goodness of all best-fittings, indicated by the coefficient of determination
$R^{2}$.

We emphasize that in Ref. [22], the hierarchical states were introduced as a
new class of eigenstates in mixed-type systems, distinct from the regular and
chaotic ones. These states, observed in the quantum kicked rotor with much
simpler phase space structure, predominantly live near the regular islands,
with only a minor contribution in the main part of the chaotic sea. They also
demonstrate a power-law decay $\hbar^{\alpha}\ (\alpha>0)$, attributed to the
transport properties associated with stickiness between regular islands and
the chaotic sea. From this definition, the hierarchical states mostly are
highly localized states as shown in Fig. 5(b1), (c1), (d1) (or the same
labeled states in Fig. 10-11). While the overlap index captures all the mixed
states across different degrees of localization, it is influenced also by the
flooding of chaotic sea to regular islands or the tunneling between regular
islands, not directly related to the effects of classical stickiness. So it
would be interesting to study the behavior of mixed states in systems without
classical stickiness, like certain lemon billiards [66, 67].

##  V Conclusions

In this work, we employed the Krylov subspace technique to generate spin
coherent states, thereby extending the investigation of the statistics of
Husimi functions of eigenstates in the kicked top, a prototype model for
studying quantum chaos. Notably, we explored system sizes up to $j=2^{13}$, a
significant advancement compared to previous studies of Husimi functions
limited to $j$ in the range of a few hundred. We examined the transition to
chaos in both classical and quantum case. All the indicators show a good
quantum-classical correspondence, specifically the mean Wehrl entropy
localization measure in the statistics of eigenstates. With this new spin
coherent states generating method, we also demonstrated that the mean Wehrl
entropy localization measure of the fully chaotic quantum kicked top, which
follows the COE spectral statistics, approaches the CUE value as the system
size $j$ increases, which is surprising [50].

More importantly, in the case of a mixed-type kicked top where regular islands
and chaotic sea coexist in the classical compact phase space, we examined
mixed eigenstates using a joint distribution of the overlap index $M$ and
Wehrl localization measure $\mathcal{L}$. $M=+1$ indicates chaotic states,
while $M=-1$ indicates regular ones. We extended the study of mixed
eigenstates with $|M|<1$ in our previous investigation [26] from $j=400$ to
$j=2^{13}$, confirming the power-law decay of the fraction of mixed states
across two orders of magnitude in $j$. This provides supporting evidence for
the PUSC and the Berry-Robnik picture in the semiclassical limit. Based on
this result from quantum kicked top, and similar results from other models
[23, 58, 65], we conjecture that the power-law decay is a universal property
for mixed eigenstates. This hypothesis merits further theoretical scrutiny.

##  VI Acknowledgement

This work was supported by the Slovenian Research and Innovation Agency (ARIS)
under the grant J1-4387.

##  Appendix A Mean ELM in integrable QKT

In Fig. 3 we have shown that the mean Wehrl ELM $\langle\mathcal{L}\rangle$ of
KT maintains a value equal to the mean Wehrl ELM of the unperturbed
Hamiltonian $\langle\mathcal{L}\rangle_{J_{x}}$, for $\gamma\leq 2$, which is
almost integrable. Fig. 9 shows that at $\gamma=1$, the mean Wehrl entropy
$\langle S\rangle\sim\frac{1}{2}\ln N$ and accordingly in the inset
$\langle\mathcal{L}\rangle\sim N^{-1/2}$. When contrasting with the uniformly
chaotic systems where $\langle S\rangle\sim\ln N$ and
$\langle\mathcal{L}\rangle\sim 1$, the distinction between integrable and
chaotic cases becomes apparent (see Fig. 4 for the behavior of $\langle
S\rangle$ and $\langle\mathcal{L}\rangle$ at $\gamma=8$). In integrable
regimes, eigenstates tend to cover a thin band around a circle on the sphere,
whereas in chaotic regimes, eigenstates exhibit delocalization across the
entire sphere [24].

Figure 9: Mean Wehrl entropy $\langle S\rangle$, versus $N=2j+1$ for QKT with
$\gamma=1$ and $\alpha=11\pi/19$. The inset is for mean Wehrl ELM
$\langle\mathcal{L}\rangle$. The dashed lines indicate the best fittings.

To derive the scaling of the mean Wehrl entropy of $J_{x}$, it is first
necessary to prove that $\langle S\rangle_{J_{x}}=\langle S\rangle_{J_{z}}$.
While eigenvectors $|j,m\rangle_{x}$ of $J_{x}$ is a unitary transformation of
eigenvectors of $J_{z}$, i.e. the Dicke basis $|j,m\rangle$, as

$\displaystyle|j,m\rangle_{x}$ $\displaystyle=\exp(-i\pi J_{y}/2)|j,m\rangle$
$\displaystyle=R(-\pi/2,0)|j,m\rangle=\sum_{m^{\prime}}d_{m^{\prime}m}^{j}(\pi/2)|m^{\prime}\rangle,$
(21)

where $R(\theta,\phi)$ is given in Eq. (11a), $d_{mm^{\prime}}^{j}(\theta)$ is
the Wigner $d$-matrix. Moreover, $R(\pi,0)|j,m\rangle=|j,-m\rangle$ while
$R(\pi,0)K$ is the conventional time-reversal operator and $K$ is the complex
conjugation on Dicke basis. The summation of the mean $q$-moment of Husimi
functions of all the eigenvectors of $J_{x}$ can then be written as

$\displaystyle\frac{1}{N}\sum_{m}W_{|m\rangle_{x}}^{(q)}=\frac{1}{4\pi}\sum_{m}\int
d\Omega|\langle\theta,\phi|j,m\rangle_{x}|^{2q}$
$\displaystyle=\frac{1}{4\pi}\sum_{m}\int d\Omega\left|\langle
j,j|R^{\dagger}(\theta,\phi)R(-\frac{\pi}{2},0)|j,m\rangle\right|^{2q},$ (22)

where $d\Omega=\sin\theta d\theta d\phi$, and $N=2j+1$. Following Eq. (11b),
we can further expand Eq. (A) with the properties of the Wigner $d$-matrix:
$d_{mm^{\prime}}^{j}(\theta)=d_{m^{\prime}m}^{j}(-\theta)=(-1)^{m^{\prime}-m}d_{m^{\prime}m}^{j}(\theta)$,
for integer $q$. It yields

$\displaystyle\frac{1}{4\pi}\sum_{m}\int
d\Omega\big{|}(1+|\xi|^{2})^{-j}\sum_{m^{\prime}}\binom{2j}{j+m^{\prime}}^{\frac{1}{2}}\xi^{j-m^{\prime}}d_{m^{\prime}m}^{j}(\frac{\pi}{2})\big{|}^{2q}$
$\displaystyle=\frac{1}{4\pi}\sum_{m}\int
d\Omega\big{|}(1+|\xi|^{2})^{-j}\sum_{m^{\prime}}\binom{2j}{j+m^{\prime}}^{\frac{1}{2}}\xi^{j-m^{\prime}}\big{|}^{2q}$
$\displaystyle=\frac{1}{4\pi}\sum_{m}\int
d\Omega|\langle\theta,\phi|j,m\rangle|^{2q}=\frac{1}{N}\sum_{m}W_{|m\rangle}^{(q)}.$
(23)

Considering the equivalence of mean $q$-moments, we establish that $\langle
S\rangle_{J_{x}}=\langle S\rangle_{J_{z}}$, as
$S^{(q)}_{|\psi\rangle}=\frac{1}{1-q}\ln W^{(q)}_{|\psi\rangle}$. While for
Dicke states, it was proved that [68]

$\displaystyle S_{|m\rangle}=\frac{2j}{2j+1}-\ln\binom{2j}{j-m}+2j\Psi(2j+1)$
$\displaystyle-(j+m)\Psi(j+m+1)-(j-m)\Psi(j-m+1),$ (24)

where $\Psi(x)$ denotes the digamma function, and obviously
$S_{|m\rangle}=S_{|-m\rangle}$. Therefore,

$\displaystyle\langle$ $\displaystyle
S\rangle_{J_{z}}=\frac{1}{N}\sum_{m=-j}^{j}S_{|m\rangle}=\frac{2j}{2j+1}-\ln(2j)!+2jB_{2j}$
$\displaystyle+\frac{2}{2j+1}\sum_{m}\left[\ln(j+m)!-(j+m)B_{j+m}\right],$
(25)

where $B_{n}=\sum_{k=1}^{n}\frac{1}{k}\approx\ln n+\gamma_{e}$. Applying
Stirling’s approximation to factorials $\ln n!\approx n\ln n+\frac{1}{2}n-n$,
for $j\to\infty$, one then gets $\langle S\rangle_{J_{z}}\sim\frac{1}{2}\ln
N$. Consequently, the mean ELM $\langle\mathcal{L}\rangle$ scales as
$N^{-1/2}$, the same scaling as the width of effective Planck cell
$j^{-1/2}\propto N^{-1/2}$.

##  Appendix B Gallery of states

Here, we provide additional examples of the joint distribution
$(\mathcal{L},M)$ for mixed-type QKT, along with representative Husimi
functions of eigenstates from different regions of the parameter space, at two
kicking strengths $\gamma=2.3$ and $3$, accompanied by corresponding SALI
plots.

Figure 10: Analogous data as in Fig. 5, but for $\gamma=2.3$. The darkest
blues in the plots of Husimi functions from (a1)-(d3) show the area where
$Q_{n}<10^{-16}$. The most localized state (d3) precisely located at the
(trivial) fixed point $(\theta,\phi)=(\pi/2,0)$, while the state (a1) of the
maximum overlap index, is strongly localized around and avoids the outer tori
of the regular region.  Figure 11: Analogous data as in Fig. 5 and 10, but now
for $\gamma=3$. Compared with $\gamma=2.3$, in this case, the fixed point at
$(\theta,\phi)=(\pi/2,0)$ becomes unstable, which is clearly shown in the SALI
plot. Interestingly, the most localized state (d3) occupies identical fixed
points, akin to when $\gamma=2.6$.

## References

  * Percival [1973] I. Percival, Journal of Physics B: Atomic and Molecular Physics 6, L229 (1973).
  * Berry [1977] M. V. Berry, Journal of Physics A: Mathematical and General 10, 2083 (1977).
  * Lichtenberg and Lieberman [2013] A. J. Lichtenberg and M. A. Lieberman,  _Regular and chaotic dynamics_ , Vol. 38 (Springer Science & Business Media, 2013).
  * Robnik [1998] M. Robnik, Nonlinear Phenomena in Complex Systems 1, 1 (1998).
  * Srednicki [1994] M. Srednicki, Physical Review E 50, 888 (1994).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016).
  * Wang  _et al._ [2022] J. Wang, M. H. Lamann, J. Richter, R. Steinigeweg, A. Dymarsky, and J. Gemmer, Physical Review Letters 128, 180601 (2022).
  * Berry and Tabor [1977] M. V. Berry and M. Tabor, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356, 375 (1977).
  * Casati  _et al._ [1980] G. Casati, F. Valz-Gris, and I. Guarnieri, Lettere al Nuovo Cimento (1971-1985) 28, 279 (1980).
  * Bohigas  _et al._ [1984] O. Bohigas, M.-J. Giannoni, and C. Schmit, Physical review letters 52, 1 (1984).
  * Müller  _et al._ [2004] S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Physical review letters 93, 014103 (2004).
  * Berry and Robnik [1984] M. V. Berry and M. Robnik, Journal of Physics A: Mathematical and General 17, 2413 (1984).
  * Tomsovic and Ullmo [1994] S. Tomsovic and D. Ullmo, Physical Review E 50, 145 (1994).
  * Bohigas  _et al._ [1993] O. Bohigas, S. Tomsovic, and D. Ullmo, Physics Reports 223, 43 (1993).
  * Doron and Frischat [1995] E. Doron and S. D. Frischat, Physical review letters 75, 3661 (1995).
  * Bäcker  _et al._ [2005] A. Bäcker, R. Ketzmerick, and A. G. Monastra, Physical review letters 94, 054102 (2005).
  * Brodier  _et al._ [2001] O. Brodier, P. Schlagheck, and D. Ullmo, Physical Review Letters 87, 064101 (2001).
  * Bäcker  _et al._ [2008a] A. Bäcker, R. Ketzmerick, S. Löck, and L. Schilling, Physical review letters 100, 104101 (2008a).
  * Bäcker  _et al._ [2008b] A. Bäcker, R. Ketzmerick, S. Löck, M. Robnik, G. Vidmar, R. Höhmann, U. Kuhl, and H.-J. Stöckmann, Physical review letters 100, 174103 (2008b).
  * Löck  _et al._ [2010] S. Löck, A. Bäcker, R. Ketzmerick, and P. Schlagheck, Physical review letters 104, 114101 (2010).
  * Ketzmerick  _et al._ [2000] R. Ketzmerick, L. Hufnagel, F. Steinbach, and M. Weiss, Physical Review Letters 85, 1214 (2000).
  * Lozej  _et al._ [2022] Č. Lozej, D. Lukman, and M. Robnik, Physical Review E 106, 054203 (2022).
  * Haake [2013] F. Haake,  _Quantum Signatures of Chaos_ , Vol. 54 (Springer Science & Business Media, 2013).
  * Haake  _et al._ [1987] F. Haake, M. Kuś, and R. Scharf, Zeitschrift für Physik B Condensed Matter 65, 381 (1987).
  * Wang and Robnik [2023] Q. Wang and M. Robnik, Physical Review E 108, 054217 (2023).
  * Sidje [1998] R. B. Sidje, ACM Transactions on Mathematical Software (TOMS) 24, 130 (1998).
  * Moler and Van Loan [2003] C. Moler and C. Van Loan, SIAM review 45, 3 (2003).
  * Golub and Van Loan [2013] G. H. Golub and C. F. Van Loan,  _Matrix computations_ (JHU press, 2013).
  * Saad [2003] Y. Saad,  _Iterative methods for sparse linear systems_ (SIAM, 2003).
  * Skokos [2001] C. Skokos, Journal of Physics A: Mathematical and General 34, 10029 (2001).
  * Skokos  _et al._ [2003] C. Skokos, C. Antonopoulos, T. C. Bountis, and M. N. Vrahatis, Progress of Theoretical Physics Supplement 150, 439 (2003).
  * Skokos  _et al._ [2004] C. Skokos, C. Antonopoulos, T. Bountis, and M. Vrahatis, Journal of Physics A: Mathematical and General 37, 6269 (2004).
  * Bountis and Skokos [2012] T. Bountis and H. Skokos,  _Complex hamiltonian dynamics_ , Vol. 10 (Springer Science & Business Media, 2012).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Physical Review B 75, 155111 (2007).
  * Atas  _et al._ [2013] Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Physical review letters 110, 084101 (2013).
  * D’Alessio and Rigol [2014] L. D’Alessio and M. Rigol, Physical Review X 4, 041048 (2014).
  * Wehrl [1978] A. Wehrl, Reviews of Modern Physics 50, 221 (1978).
  * Perelomov [1972] A. M. Perelomov, Communications in Mathematical Physics 26, 222 (1972).
  * Klauder and Skagerstam [1985] J. R. Klauder and B.-S. Skagerstam,  _Coherent states: applications in physics and mathematical physics_ (World scientific, 1985).
  * Zhang  _et al._ [1990] W.-M. Zhang, R. Gilmore, _et al._ , Reviews of Modern Physics 62, 867 (1990).
  * Wehrl [1979] A. Wehrl, Reports on Mathematical Physics 16, 353 (1979).
  * Lieb [1978] E. H. Lieb, Communications in Mathematical Physics 62, 35 (1978).
  * Lieb and Solovej [1991] E. H. Lieb and J. P. Solovej, letters in mathematical physics 22, 145 (1991).
  * Lieb and Solovej [2014] E. H. Lieb and J. P. Solovej, Acta Mathematica 212, 379 (2014).
  * Lieb and Solovej [2016] E. H. Lieb and J. P. Solovej, Communications in Mathematical Physics 348, 567 (2016).
  * Schupp [2022] P. Schupp, arXiv preprint arXiv:2203.08095 (2022).
  * Wootters [1990] W. K. Wootters, Foundations of Physics 20, 1365 (1990).
  * Jones [1990] K. R. Jones, Journal of Physics A: Mathematical and General 23, L1247 (1990).
  * Gnutzmann and Zyczkowski [2001] S. Gnutzmann and K. Zyczkowski, Journal of Physics A: Mathematical and General 34, 10123 (2001).
  * Bogomolny [1992] E. B. Bogomolny, Nonlinearity 5, 805 (1992).
  * Prosen [1996] T. Prosen, Physica D: Nonlinear Phenomena 91, 244 (1996).
  * Batistić and Robnik [2013a] B. Batistić and M. Robnik, Journal of Physics A: Mathematical and Theoretical 46, 315102 (2013a).
  * Batistić  _et al._ [2013] B. Batistić, T. Manos, and M. Robnik, Europhysics Letters 102, 50008 (2013).
  * Batistić and Robnik [2013b] B. Batistić and M. Robnik, Physical Review E 88, 052913 (2013b).
  * Batistić  _et al._ [2019] B. Batistić, Č. Lozej, and M. Robnik, Physical Review E 100, 062208 (2019).
  * Batistić  _et al._ [2020] B. Batistić, Č. Lozej, and M. Robnik, Nonlinear Phenomena in Complex Systems 23, 17 (2020).
  * Yan and Robnik [2024] H. Yan and M. Robnik, arXiv preprint arXiv:2401.13070 (2024).
  * MacKay  _et al._ [1984] R. MacKay, J. Meiss, and I. Percival, Physica D: Nonlinear Phenomena 13, 55 (1984).
  * Geisel  _et al._ [1986] T. Geisel, G. Radons, and J. Rubner, Physical review letters 57, 2883 (1986).
  * Bohigas  _et al._ [1990] O. Bohigas, S. Tomsovic, and D. Ullmo, Physical review letters 64, 1479 (1990).
  * Michler  _et al._ [2012] M. Michler, A. Bäcker, R. Ketzmerick, H.-J. Stöckmann, and S. Tomsovic, Physical Review Letters 109, 234101 (2012).
  * Stöber  _et al._ [2024] J. Stöber, A. Bäcker, and R. Ketzmerick, Physical Review Letters 132, 047201 (2024).
  * Meiss [2015] J. Meiss, Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (2015).
  * Wang and Robnik [2024] Q. Wang and M. Robnik, Physical Review E 109, 024225 (2024).
  * Lozej  _et al._ [2021a] C. Lozej, D. Lukman, and M. Robnik, Nonlinear Phenomena in Complex Systems 24, 1 (2021a).
  * Lozej  _et al._ [2021b] Č. Lozej, D. Lukman, and M. Robnik, Physical Review E 103, 012204 (2021b).
  * Lee [1988] C. T. Lee, Journal of Physics A: Mathematical and General 21, 3749 (1988).


# NAAS: Neural Accelerator Architecture Search

Yujun Lin* MIT Mengtian Yang* SJTU
https://tinyml.mit.edu Song Han *Equally contributed to this work. MIT

###### Abstract

Data-driven, automatic design space exploration of neural accelerator
architecture is desirable for specialization and productivity. Previous
frameworks focus on sizing the numerical architectural hyper-parameters while
neglect searching the PE connectivities and compiler mappings. To tackle this
challenge, we propose Neural Accelerator Architecture Search (NAAS) that
holistically searches the neural network architecture, accelerator
architecture and compiler mapping in one optimization loop. NAAS composes
highly matched architectures together with efficient mapping. As a data-driven
approach, NAAS rivals the human design Eyeriss by 4.4$\times$ EDP reduction
with 2.7% accuracy improvement on ImageNet under the same computation
resource, and offers 1.4$\times$ to 3.5$\times$ EDP reduction than only sizing
the architectural hyper-parameters.

††publicationid: pubid: 978-1-6654-3274-0/21/$31.00 ©2021 IEEE

##  I Introduction

Neural architecture and accelerator architecture co-design is important to
enable specialization and acceleration. It covers three aspects: designing the
neural network, designing the accelerator, and the compiler that maps the
model on the accelerator. The design space of each dimension is listed in
Table I, with more than $10^{800}$ choices for a 50-layer neural network.
Given the huge design space, data-driven approach is desirable, where new
architecture design evolves as new designs and rewards are collected. Recent
work on hardware-aware neural architecture search (NAS) and auto compiler
optimization have successfully leverage machine learning algorithms to
automatically explore the design space. However, these works only focuses on
off-the-shelf hardware [1, 2, 3, 4, 5, 6], and neglect the freedom in the
hardware design space.

The interactions between the neural architecture and the accelerator
architecture is illustrated in Table II. The correlations are complicated and
vary from hardware to hardware: for instance, tiled input channels should be
multiples of #rows of compute array in NVDLA, while #rows is related to the
kernel size in Eyeriss. It is important to consider all the correlations and
make them fit. A tuple of perfectly matched neural architecture, accelerator
architecture, and mapping strategy will improve the utilization of the compute
array and on-chip memory, maximizing efficiency and performance.

The potential of exploring both neural and accelerator architecture has been
proven on FPGA platforms [7, 8, 9, 10] where HLS is applied to generate FPGA
accelerator. Earlier work on accelerator architecture search [11, 12, 13] only
search the architectural sizing while neglecting the PE connectivity (_e.g_.,
array shape and parallel dimensions) and compiler mappings, which impact the
hardware efficiency.

TABLE I: Neural-Accelerator architecture search space.

Accelerator | Compute Array Size ($\\#$rows/$\\#$columns)
---|---
(Input/Weight/Output) Buffer Size
PE Inter-connection (Dataflow)
Compiler Mapping | Loop Order, Loop Tiling Sizes
Neural Network | #Layers, #Channels, Kernel Size
Block Structure, Input Size

TABLE II: The complicated correlation between neural and accelerator design
space. It differs from accelerator to accelerator: N is NVDLA and E is
Eyeriss.

Accelerator | Neural Architecture Design Space
---|---
Parameter | Input | Output | Kernel | Feature
Space | Channels | Channels | Size | Map Size
Array #rows | N |  | E |
Array #cols |  | N |  | E
IBUF Size | N/E |  |  | N/E
WBUF Size | N/E | N/E | N/E |
OBUF Size |  | N/E |  | N/E

We push beyond searching only hardware hyper-parameters and propose the Neural
Accelerator Architecture Search (NAAS), which fully exploits the hardware
design space and compiler mapping strategies at the same time. Unlike prior
work [11] which formulate the hardware parameter search as a pure sizing
optimization, NAAS models the co-search as a two-level optimization problem,
where each level is a combination of indexing, ordering and sizing
optimization. To tackle such challenges, we propose an encoding method which
is able to encode the non-numerical parameters such as loop order and parallel
dimension chosen as numerical parameters for optimization. As shown in Figure
1, the outer loop of NAAS optimizes the accelerator architecture while the
inner loop optimizes the compiler mapping strategies.

Combining both spaces greatly enlarges the optimization space: within the same
#PEs and on-chip memory resources as EdgeTPU there are at least $10^{11}$
hardware candidates and $10^{17}$ mapping candidates for each layer, which
composes $10^{11+50\cdot 17}=10^{861}$ possible combinations in the joint
search space for ResNet-50, while there are only $10^{4}$ hardware candidates
in NASAIC’s design space. To efficiently search over the large design space,
NAAS leverages the biologically-inspired evolution-based algorithm rather than
meta-controller-based algorithm to improve the sample efficiency. It keeps
improving the quality of the candidate population by ruling out the inferior
and generating from the fittest. Thanks to the low search cost, NAAS can be
easily integrated with hardware-aware NAS algorithm by adding another
optimization level (Figure 1), achieving the joint search.

Extensive experiments verify the effectiveness of our framework. Under the
same #PE and on-chip memory constraints, the NAAS is able to deliver
2.6$\times$, 4.4$\times$ speedup and 2.1$\times$, 1.4$\times$ energy savings
on average compared to Eyeriss [14], NVDLA [15] design respectively.
Integrated with Once-For-All NAS algorithm [4], NAAS further improves the
top-1 accuracy on ImageNet by 2.7% without hurting the hardware performance.
Using the similar compute resources, NAAS achieves 3.0$\times$, 1.9$\times$
EDP improvements compared to Neural-Hardware Architecture Search [12], and
NASAIC [11] respectively.

##  II Neural Accelerator Architecture Search

Figure 1: Neural Accelerator Architecture Search.

Figure 1 shows the optimization flow of Neural Accelerator Architecture Search
(NAAS). NAAS explores the design space of accelerators, and compiler’s
mappings simultaneously.

###  II-A Accelerator Architecture Search

Figure 2: Encoding the accelerator design and compiler mapping into vectors.
An accelerator design is described by architectural sizing and connectivity
parameters, where PE inter-connection is presented by parallel dimension. A
mapping strategy is represented by loop order and corresponding tiling at each
dimension of array (the innermost level is PE level). The mapping description
in the figure is in MAESTRO[16] format which fuses the array parameters and
mapping strategy.  Figure 3: We convert the non-numerical ordering
optimization into numerical optimization for parallel dimensions and loop
order searching. Both accelerator optimizer and mapping optimizer will
generate an importance value for each dimension. The left part shows the
parallel dimensions of this 2D array candidate are the dimensions with the
largest two importance value. The right part shows the order of each dimension
in for-loop is in the decreasing order of its importance value. Figure 4: The
average of EDP decreases as NAAS is learning.

#### Design Space

The accelerator design knobs can be categorized into two classes:

  1. 1.

Architectural Sizing: the number of processing elements (#PEs), private
scratch pad size (L1 size), global buffer size (L2 size), and memory
bandwidth.

  2. 2.

Connectivity Parameters: the number of array dimensions (1D, 2D or 3D array),
array size at each dimension, and the inter-PE connections.

Most state-of-art searching frameworks only contains architectural sizing
parameters in their design space. These sizing parameters are numerical and
can be easily embedded into vectors during search. On the other hand, PE
connectivity is difficult to encode as vectors since they are not numerical
numbers. Moreover, changing the connectivity requires re-designing the
compiler mapping strategies, which extremely increase the searching cost. In
NAAS, besides the architectural sizing parameters which are common in other
frameworks, we introduce the connectivity parameters into our search space,
making it possible to search among 1D, 2D and 3D array as well, and thus our
design space includes almost the entire accelerator design space for neural
network accelerators.

#### Encoding

We first model the PE connectivity as the choices of parallel dimensions. For
example, parallelism in input channels (C) means a reduction connection of the
partial sum register inside each PE. Parallelism in output channels means a
broadcast to input feature register inside each PE. The most straight-forward
method to encode the parallel dimension choice is to enumerate all possible
parallelism situations and choose the index of the enumeration as the encoding
value. However, since the increment or decrement of indexes does not convey
any physical information, it is hard to be optimized.

To solve this problem, we proposed the “importance-based” encoding method for
choosing parallelism dimensions in the dataflow and convert the indexing
optimization into the sizing optimization. For each dimension, our optimizer
will generate an importance value. To get the corresponding parallel
dimensions, we first collect all the importance value, then sort them in
decreasing order, and select the first k dimensions as the parallel dimensions
of a k-D compute array. As shown in the left of Figure 3, the generated
candidate is a 2D array with size $16\times 16$. To find the parallel
dimension for this 2D array candidate, The importance values are first
generated for 6 dimensions in the same way as other numerical parameters in
the encoding vector. We then sort the value in decreasing order and determine
the new order of the dimensions. Since the importance value of “C” and “K” are
the largest two value, we finally select “C” and “K” as the parallel
dimensions of this 2D array. The importance value of the dimension represents
the priority of the parallelism: a larger value indicates a higher priority
and a higher possibility to be paralleled in the computation loop nest, which
contains higher relativity with accelerator design compared to indexes of
enumerations.

For other numerical parameters, we use the straight-forward encoding method.
The whole hardware encoding vector is shown in Figure 2, which contains all of
the necessary parameters to represent an accelerator design paradigm.

#### Evolution Search

We leverage the evolution strategy [17] to find the best solution during the
exploration. In order to take both latency and energy into consideration, we
choose the widely used metric Energy-Delay Product (EDP) to evaluate a given
accelerator configuration on a specific neural network workload. At each
evolution iteration, we first sample a set of candidates according to a
multivariate normal distribution in $\left[0,1\right]^{\left|\theta\right|}$.
Each candidate is represented as a $\left|\theta\right|$-dimension vector.
These candidates are then projected to hardware encoding vectors and decoded
into accelerator design. We rule out the invalid accelerator samples and keep
sampling until the candidate set reaches a predefined size (population size)
in our experiments. To evaluate the candidate performance, we need to perform
mapping strategy search in Section II-B on each benchmark and adopt the best
searched result as the EDP reward of this candidate. After evaluating all the
candidates on the benchmarks, we update the sampling distribution based on the
relative ordering of their EDP. Specifically, we select the top solutions as
the “parents” of the next generation and use their center to generate the new
mean of the sampling distribution. We update the covariance matrix of the
distribution to increase the likelihood of generating samples near the parents
[17]. We then repeat such sampling and updating process.

Figure 4 shows the statistics of energy-delay products of hardware candidates
in the population. As the optimization continues, the EDP mean of NAAS
candidates decreases while that of random search remains high, which indicates
that NAAS gradually improves the range of hardware selections.

###  II-B Compiler Mapping Strategy Search

The performance and energy efficiency of deep learning accelerators also
depend on how to map the neural network task on the accelerator. The search
space of compiler mapping strategy is much larger than accelerator design,
since different convolution layers may not share the same optimal mapping
strategy. Hence we optimize the mapping for each layer independently using the
similar evolution-based search algorithm to accelerator design search in
Section II-A.

The compiler mapping strategy consists of two components: the execution order
and the tiling size of each for-loop dimension. Similar to the accelerator
design search, the order of for-loop dimensions is non-trivial. Rather than
enumerating all of the possible execution orders and using indexes as
encoding, we use the similar “importance-based” encoding methods in Section
II-A. For each dimension of the array (corresponding to each level of for-loop
nests), the mapping optimizer will assign each convolution dimension with an
importance value. The dimension with the highest importance will become the
outermost loop while the one with the lowest importance will be placed at the
innermost in the loop nests. The right of Figure 3 gives an example. The
optimizer firstly generates the importance values for 6 dimensions, then sort
the value in decreasing order and determine the corresponding order of the
dimensions. Since “C” and “R” dimension have largest value 5, they will become
the outermost loops. “S” dimension has the smallest value 1, so it is the
innermost dimension in the loops. This strategy is interpretable, since the
importance value represents the data locality of the dimension: the dimension
labeled as most important has the best data locality since it is the outermost
loop, while the dimension labeled as least important has the poorest data
locality therefore it is the innermost loop.

As for tiling sizes, since they are highly related to the network parameters,
we use the scaling ratio rather than the absolute tiling value. Hence, the
tiling sizes are still numerical parameters and able to adapt to different
networks. The right part of figure 3 illustrates the composition of the
mapping encoding vector. Note that for PE level we need to ensure that there
is only one MAC in a PE, so we only search the loop order at PE level. For
each array level, the encoding vector contains both the execution order and
tiling size for each for-loop dimension.

###  II-C Integrated with Neural Architecture Search

Thanks to the low search cost, we can integrate our framework with neural
architecture search to achieve neural-accelerator-compiler co-design. Figure 1
illustrates integrating NAAS with NAS. The joint design space is huge, and in
order to improve the search efficiency, we choose and adapt Once-For-All NAS
algorithm for NAAS. First, NAAS generates a pool of accelerator candidates.
For each accelerator candidate, we sample a network candidate from NAS
framework which satisfies the pre-defined accuracy requirement. Since each
subnet of Once-For-All network is well trained, the accuracy evaluation is
fast. We then apply the compiler mapping strategy search for the network
candidate on the corresponding accelerator candidate. NAS optimizer will
update using the searched EDP as a reward. Until NAS optimizer reaches its
iteration limitations, and feedback the EDP of the best network candidate to
accelerator design optimizer. We repeated the process until the best-fitted
design is found. In the end, we obtain a tuple of matched accelerator, neural
network, and its mapping strategy with guaranteed accuracy and lowest EDP.

##  III Evaluation

We evaluate Neural Accelerator Architecture Search’s performance improvement
step by step: 1) the improvement from applying NAAS given the same hardware
resource; 2) performance of NAAS which integrates the Once-For-All to achieve
the neural-accelerator-compiler co-design.

Figure 5: Speedup and energy savings when NAAS searches the accelerator and
mapping for a set of network benchmarks within the baseline hardware
resources. For instance, using the same #PEs and SRAM size as EdgeTPU, the
accelerator designed by NAAS runs VGG16, ResNet50 and UNet 2.6$\times$ faster
than EdgeTPU on average. Figure 6: Speedup and energy savings when NAAS
searches the accelerator and mapping for single network benchmark within the
baseline hardware resources. Figure 7: NAAS offers different architectures for
different networks when given different computation resources: (a) 2D array
with K-X’ parallelism for ResNet using Eyeriss resource; (b) 2D array with
C-X’ parallelism for VGG16 using EdgeTPU resource; (c) 3D array with C-K-X’
parallelism for VGG16 using Shidiannao resource. Figure 8: Compared to
searching the architectural sizing only [11, 12], searching the connectivity
parameters and mapping strategies as well achieves considerable EDP reduction.
Figure 9: Compared to index-based encoding, importance-based encoding achieves
better EDP reduction. Figure 10: Accuracy _vs_. Normalized EDP on ImageNet
(batch = 1) under Eyeriss hardware resource. Integrating NAAS with OFA NAS
further improves model accuracy.

###  III-A Evaluation Environment

#### Design Space of NAAS

We select four different resource constraints based on EdgeTPU, NVDLA [15],
Eyeriss [14] and Shidiannao [18]. When comparing to each baseline
architecture, NAAS is conducted within corresponding computation resource
constraint including the maximum #PEs, the maximum total on-chip memory size,
and the NoC bandwidth. NAAS searches #PEs at stride of 8, buffer sizes at
stride of 16B, array sizes at stride of 2.

#### CNN Benchmarks

We select 6 widely-used CNN models as our benchmarks. The benchmarks are
divided into two sets: classic large-scale networks (VGG16, ResNet50, UNet)
and light-weight efficient mobile networks (MobileNetV2, SqueezeNet, MNasNet).
Five deployment scenarios are divided accordingly: we conduct NAAS for large
models with more hardware resources (EdgeTPU, NVDLA with 1024 PEs), and for
small models with limited hardware resources (ShiDianNao, Eyeriss, and NVDLA
with 256 PEs).

#### Design Space in Once-For-All NAS

When integrating with Once-For-All NAS, the neural architecture space is
modified from ResNet-50 design space following the open-sourced library [4].
There are 3 width multiplier choices (0.65, 0.8, 1.0) and 18 residual blocks
at maximum, where each block consists of three convolutions and has 3
reduction ratios (0.2, 0.25, 0.35). Input image size ranges from 128 to 256 at
stride of 16. In total there are $10^{13}$ possible neural architectures.

###  III-B Improvement from NAAS

Figure 5 shows the speedup and energy savings of NAAS using the same hardware
resources compared to the baseline architectures. When running large-scale
models, NAAS delivers 2.6$\times$, 2.2$\times$ speedup and 1.1$\times$,
1.1$\times$ energy savings on average compared to EdgeTPU and NVDLA-1024 .
Though NAAS tries to provide a balanced performance on all benchmarks by using
geomean EDP as reward, VGG16 workload sees the highest gains from NAAS. When
inferencing light-weight models, NAAS achieves 4.4$\times$, 1.7$\times$,
4.4$\times$ speedup and 2.1$\times$, 1.4$\times$, 4.9$\times$ energy savings
on average compared to Eyeriss, NVDLA-256, and ShiDianNao. Similar to
searching for large models, different models obtain different benefits from
NAAS under different resource constraints.

Figure 7 demonstrates three examples of searched architectures. When given
different computation resources, for different NN models, NAAS provides
different solutions beyond numerical design parameter tuning. Different
dataflow parallelisms determine the different PE micro-architecture and thus
PE connectivies and even feature/weight/partial-sum buffer placement.

Figure 8 illustrates the benefits of searching connectivity parameters and
mapping strategy compared to searching architectural sizing only. NAAS
outperforms architectural sizing search by 3.52$\times$, 1.42$\times$EDP
reduction on VGG and MobileNetV2 within EdgeTPU resources, as well as
2.61$\times$, 1.62$\times$improvement under NVDLA-1024 resources.

Figure 9 further shows the EDP reduction using different encoding methods for
non-numerical parameters in hardware and mapping encoding vectors. Our
proposed importance-based encoding method significantly improves the
performance of optimization by reducing EDP from 1.4$\times$to 7.4$\times$.

###  III-C More Improvement from Integrating NAS

Different from accelerator architectures, neural architectures have much more
knobs to tune (_e.g_., network depths, channel numbers, input image sizes),
providing us with more room to optimize. To illustrate the benefit of NAAS
with NAS, we evaluate on ResNet50 with hardware resources similar to Eyeriss.
Figure 10 shows that NAAS (accelerator only) outperforms Neural-Hardware
Architecture Search (NHAS) [12] (which only searches the neural architecture
and the accelerator architectural sizing) by 3.01$\times$ EDP improvement. By
integrating with neural architecture search, NAAS achieves 4.88$\times$ EDP
improvement in total as well as 2.7% top-1 accuracy improvement on ImageNet
dataset than Eyeriss running ResNet50.

#### Comparison to NASAIC

We also compare our NAAS with previous work NASAIC [11] in Table III. NASAIC
adopts DLA [15], ShiDianNao [18] as source architectures for its heterogeneous
accelerator, and only searches the allocation of #PEs and NoC bandwidth. In
contrast, we explore more possible improvements by adapting both accelerator
design searching and mapping strategy searching. Inferencing the same network
searched by NASAIC, NAAS outperforms NASAIC by 1.88$\times$ EDP improvement
(3.75$\times$ latency improvement with double energy cost) using the same
design constraints.

TABLE III: NAAS (accelerator only) achieves better EDP.

Search | Arch | Cifar-10 | Latency | Energy | EDP
---|---|---|---|---|---
Approach | Accuracy | (cycles) | (nJ) | (cycles$\cdot$nJ)
NASAIC [11] | DLA [15] | 93.2 | 3e5 | 1e9 | 3e14
Shi [18] | 91.1
NAAS | DLA [15] | 93.2 | 8e4 | 2e9 | 2e14

TABLE IV: We achieve much lower search cost on ImageNet (Gds: GPU days. N:
number of development scenarios.)

Approach | Co-Search | NN Training | Total | AWS | CO2
---|---|---|---|---|---
Cost (Gds) | Cost (Gds) | Cost (Gds) | Cost | Emission
NASAIC | 500$\times$12N = 6000N | 16N | 6000N | $441,000N | 41,000N lbs
NHAS | 12 + 4N | 16N | 12 + 20N | $1,500N | 150N lbs
Ours |  $<$ 0.25N | 50 |  $<$ 50 + 0.25N |  $<$ $18N |  $<$ 2N lbs
* NASAIC’s search cost is an optimistic projection from Cifar10.
* AWS cost $75/Gd and CO2 Emission is 7.5 lbs/Gd.

#### Search Cost

Table IV reports the search cost of our NAAS compared to NASAIC and NHAS when
developing accelerator and network for $N$ development scenarios, where “AWS
Cost” is calculated based on the price of on-demand P3.16xlarge instances, and
“CO2 Emission” is calculated based on Strubell _et al_. [19]. NASAIC relies on
a meta-controller-based search algorithm which requires training every neural
architecture candidates from scratch. NHAS [12] decouples training and
searching in NAS but it also requires retraining the searched network for each
deployment which costs 16$N$ GPU days. The high sample efficiency of NAAS
makes it possible to integrate the Once-For-All NAS in one search loop. As a
conservative estimation, NAAS saves more than 120$\times$ search cost compared
to NASAIC on ImageNet.

##  IV Related Works

#### Accelerator Design-Space Exploration

Ealier work focuses on fine-grained hardware resource assignment for
deployment on FPGAs [7, 20, 21, 22, 10, 13]. Several work focuses on co-
designing neural architectures and ASIC accelerators. Yang  _et al_. [11]
(NASAIC) devise a controller that simultaneously predicts neural architectures
as well as the selection policy of various IPs in a heterogeneous accelerator.
Lin  _et al_. [12] focuses on optimizing the micro architecture parameters
such as the array size and buffer size for given accelerator design while
searching the quantized neural architecture. Besides, some work focuses on
optimizing compiler mapping strategy on fixed architectures. Chen  _et al_.
[23] (TVM) designed a searching framework for optimizing for-loop execution on
CPU/GPU or other fixed platforms. Mu  _et al_. [24] proposed a new searching
algorithm for tuning mapping strategy on fixed GPU architecture. On the
contrary, our work explores not only the sizing parameters but also the
connectivity parameters and the compiler mapping strategy. We also explore the
neural architecture space to further improve the performance. Plenty of work
provides the modeling platform for design space exploration [25, 26, 27, 28].
We choose MAESTRO [28] as the accelerator evaluation backend.

#### AutoML and Hardware-Aware NAS

Researchers have looked to automate the neural network design using AutoML.
Grid search [28, 29] and reinforcement learning with meta-controller [30, 2,
11] both suffer from prohibitive search cost. One-shot-network-based
frameworks [1, 3, 31] achieved high performance at a relatively low search
cost. These NAS algorithms require retraining the searched networks while Cai
_et al_. [4] proposed Once-For-All network of which the subnets are well
trained and can be directly extracted for deployment. Recent neural
architecture search (NAS) frameworks started to incorporate the hardware into
the search feedback loop [2, 9, 1, 8, 3, 4], though they have not explored
hardware accelerator optimization.

##  V Conclusion

We propose an evolution-based accelerator-compiler co-search framework, NAAS.
It not only searches the architecture sizing parameters but also the PE
connectivity and compiler mapping strategy. Integrated with the Once-for-All
NAS algorithm, it explores the search spaces of neural architectures,
accelerator architectures, and mapping strategies together while reducing the
search cost by 120$\times$ compared with previous work. Extensive experiments
verify the effectiveness of NAAS. Within the same computation resources as
Eyeriss [14], NAAS provides 4.4$\times$ energy-delay-product reduction with
2.7% top-1 accuracy improvement on ImageNet dataset compared to directly
running ResNet-50 on Eyeriss. Using the similar computation resources, NAAS
integrated with NAS achieves 3.0$\times$, 1.9$\times$ EDP improvements
compared to NHAS [12], and NASAIC [11] respectively.

Acknowledgements. This work was supported by NSF CAREER Award #1943349 and SRC
GRC program under task 2944.001. We also thank AWS Machine Learning Research
Awards for the computational resource.

## References

  * [1] H. Cai, L. Zhu, and S. Han, “ProxylessNAS: Direct neural architecture search on target task and hardware,” in _ICLR_ , 2019.
  * [2] Q. Lu _et al._ , “On neural architecture search for resource-constrained hardware platforms,” 2019.
  * [3] K. Wang _et al._ , “HAQ: Hardware-aware automated quantization with mixed precision,” in _CVPR_ , 2019, pp. 8612–8620.
  * [4] H. Cai _et al._ , “Once for all: Train one network and specialize it for efficient deployment,” in _ICLR_ , 2020.
  * [5] T. Chen _et al._ , “Learning to optimize tensor programs,” in _NeurIPS_ , 2018, pp. 3389–3400.
  * [6] S.-C. Kao and T. Krishna, “Gamma: Mapping space exploration via genetic algorithm,” in _ICCAD’20_ , 2020.
  * [7] C. Hao _et al._ , “FPGA/DNN co-design: An efficient design methodology for 1ot intelligence on the edge,” in _DAC_ , 2019, pp. 1–6.
  * [8] Y. Li _et al._ , “EDD: Efficient differentiable dnn architecture and implementation co-search for embedded ai solutions,” _DAC_ , 2020.
  * [9] X. Zhang _et al._ , “Skynet: a hardware-efficient method for object detection and tracking on embedded systems,” 2020.
  * [10] S.-C. Kao _et al._ , “Confuciux: Autonomous hardware resource assignment for dnn accelerators using reinforcement learning,” in _MICRO_ , 2020.
  * [11] L. Yang _et al._ , “Co-Exploration of Neural Architectures and Heterogeneous ASIC Accelerator Designs Targeting Multiple Tasks,” 2020.
  * [12] Y. Lin _et al._ , “Neural hardware architecture search,” _NeurIPS WS_ , 2019.
  * [13] W. Jiang _et al._ , “Hardware/software co-exploration of neural architectures,” _TCAD_ , 2020.
  * [14] Y.-H. Chen _et al._ , “Eyeriss: An energy-efficient reconfigurable accelerator for deep convolutional neural networks,” _IJSSC_ , vol. 52, no. 1, 2016.
  * [15] NVIDIA, “NVIDIA deep learning accelerator,” 2017. [Online]. Available: http://nvidia.org
  * [16] J. Albericio _et al._ , “Cnvlutin: Ineffectual-neuron-free deep neural network computing,” _ACM SIGARCH_ , vol. 44, no. 3, pp. 1–13, 2016.
  * [17] N. Hansen, “The CMA evolution strategy: a comparing review,” in _Towards a new evolutionary computation_. Springer, 2006, pp. 75–102.
  * [18] Z. Du _et al._ , “Shidiannao: Shifting vision processing closer to the sensor,” in _ISCA_ , 2015, pp. 92–104.
  * [19] E. Strubell _et al._ , “Energy and policy considerations for deep learning in nlp,” in _ACL_ , 2019, pp. 3645–3650.
  * [20] M. Motamedi _et al._ , “Design space exploration of fpga-based deep convolutional neural networks,” in _ASP-DAC_. IEEE, 2016, pp. 575–580.
  * [21] G. Zhong _et al._ , “Design space exploration of fpga-based accelerators with multi-level parallelism,” in _DATE_. IEEE, 2017, pp. 1141–1146.
  * [22] Y. Chen _et al._ , “Cloud-dnn: An open framework for mapping dnn models to cloud fpgas,” in _FPGA_ , 2019, pp. 73–82.
  * [23] T. Chen _et al._ , “$\\{$TVM$\\}$: An automated end-to-end optimizing compiler for deep learning,” in _13th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 18)_, 2018, pp. 578–594.
  * [24] J. Mu _et al._ , “A history-based auto-tuning framework for fast and high-performance dnn design on gpu,” in _2020 57th ACM/IEEE Design Automation Conference (DAC)_. IEEE, 2020, pp. 1–6.
  * [25] Y. S. Shao _et al._ , “The aladdin approach to accelerator design and modeling,” _IEEE Micro_ , vol. 35, no. 3, pp. 58–70, 2015.
  * [26] Y. N. Wu and V. Sze, “Accelergy: An architecture-level energy estimation methodology for accelerator designs,” in _TCAD_ , 2019.
  * [27] A. Parashar _et al._ , “Timeloop: A systematic approach to dnn accelerator evaluation,” in _ISPASS_. IEEE, 2019, pp. 304–315.
  * [28] H. Kwon _et al._ , “Understanding reuse, performance, and hardware cost of dnn dataflow: A data-centric approach,” in _ISCA_ , 2019, pp. 754–768.
  * [29] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _ICML_ , 2019, pp. 6105–6114.
  * [30] B. Zoph and Q. V. Le, “Neural architecture search with reinforcement learning,” in _ICLR_ , 2017.
  * [31] Z. Guo _et al._ , “Single path one-shot neural architecture search with uniform sampling,” _arXiv preprint_ , 2019.