††thanks: Corresponding author<EMAIL_ADDRESS>
# Supplementary Materials for “Attention to Quantum Complexity”

Hyejin Kim† Department of Physics, Cornell University, Ithaca, NY, USA  Yiqing
Zhou† Department of Physics, Cornell University, Ithaca, NY, USA  Yichen Xu†
Department of Physics, Cornell University, Ithaca, NY, USA  Kaarthik Varma
Department of Physics, Cornell University, Ithaca, NY, USA  Amir H. Karamlou
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Ilan T. Rosen  Research Laboratory of Electronics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Jesse C. Hoke  Google
Research, Mountain View, CA, USA  Department of Physics, Stanford University,
Stanford, CA, USA  Chao Wan  Department of Computer Science, Cornell
University, Ithaca, New York 14853, USA  Jin Peng Zhou  Department of Computer
Science, Cornell University, Ithaca, New York 14853, USA  William D. Oliver
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Research Laboratory of Electronics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Department of Electrical Engineering and
Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  Yuri D. Lensky  Department of Physics, Cornell University, Ithaca, NY,
USA  Google Research, Mountain View, CA, USA  Kilian Q. Weinberger  Department
of Computer Science, Cornell University, Ithaca, New York 14853, USA  Eun-Ah
Kim  Department of Physics, Cornell University, Ithaca, NY, USA  Google
Research, Mountain View, CA, USA  Department of Physics, Ewha Womans
University, Seoul, South Korea
†These authors contributed equally to this work.

###### Contents

  1. A QuAN Architecture
    1. A.1 Input
    2. A.2 Convolution
    3. A.3 Mini-set Self-attention block (MSSAB) Encoder
      1. A.3.1 Parallel SAB
      2. A.3.2 Recurrent AB (RecAB)
      3. A.3.3 Reducing AB (RedAB)
      4. A.3.4 Computational complexity of MSSAB
    4. A.4 Pooling Attention Block (PAB) Decoder
  2. B Comparison of QuAN and other ML architectures
  3. C Driven Hard-core Bose-Hubbard model
    1. C.1 Data acquisition
    2. C.2 Data preprocessing
    3. C.3 Training and testing procedure
    4. C.4 Machine learning details
  4. D Random quantum circuit
    1. D.1 Data acquisition
      1. D.1.1 Quantum processor details and experimental procedure
      2. D.1.2 Simulation and linear cross-entropy benchmarking (XEB)
    2. D.2 Data preprocessing
    3. D.3 Training and testing procedure
    4. D.4 Machine learning details
  5. E Toric code simulation
    1. E.1 Data acquisition
    2. E.2 Data preprocessing
    3. E.3 Training and testing procedure
    4. E.4 Benchmarking machine results
      1. E.4.1 Benchmarking to locally error-corrected decoration (LED)
      2. E.4.2 Benchmarking to SMLP
    5. E.5 Machine analysis: PAB as an importance-sampler

##  Appendix A QuAN Architecture

Capturing quantum fluctuation and respecting samples’ permutation invariance
is essential to learn information about quantum states from measurement
snapshots. Here, we introduce QuAN, a machine learning model which uses
attention mechanism [1] to learn inter-snapshot correlations while respecting
the permutation invariance of the snapshots [2, 3]. Learning high-order
moments can be a practical way of learning characteristics of the distribution
of bit-strings. However, direct application of existing self-attention blocks
as in Ref. [3] is inefficient and allows us to access only low-order moments.
To study complex quantum systems where high-order moments become crucial, we
propose (in Section A.3) a new encoding scheme, named mini-set self-attention
block (MSSAB), that can efficiently sample high-order correlations with low
computational cost.

In this section, we present details of the inner structure of QuAN. A detailed
structure of QuAN architecture is presented in the main text Fig. 1(e). The
rest of this section is organized as follows: First, we discuss in Section A.1
the input data structure taken by QuAN. Second, we discuss in Section A.2 the
usage of a two-dimensional convolution layer and any preprocessing of the raw
measurement snapshots before inputting into the encoder. Third, we discuss in
Section A.3 how the workhorse in QuAN - the mini-set self-attention block
(MSSAB) is used in the encoder to capture inter-snapshot correlations. Lastly,
we discuss in Section A.4 the pooling attention block (PAB) in the decoder.

###  A.1 Input

QuAN takes a set of two-dimensional bit-string of $N_{q}$ qubits as input
data, where each bit-string consists of $N_{q}$ binary values (0 or 1). We
denote the inputs to QuAN, which are $i$-th set (datapoint)
$\mathbb{X}_{i}=\\{B_{i,\alpha}\\}_{\alpha=1}^{N}$ consisting of $N$ two-
dimensional binary-valued arrays $B_{i,\alpha}$, where $\alpha$ indexes the
elements within the set and $N$ is the set size. (Each binary-valued array has
$N_{q}$ entries indexed by $\mu$, i.e. $B_{i,\alpha,\mu}=0$ or $1$. To avoid
confusion between the set element index $\alpha\in\\{1,\cdots,N\\}$ and the
spatial dimension index $\mu\in\\{1,\cdots,N_{q}\\}$, we relocate the index
$\alpha$ from superscript to subscript, s.t. $B^{\alpha}_{i,\mu}\equiv
B_{i,\alpha,\mu}$.) Once $\mathbb{X}_{i}$ is inputted into QuAN, the output is
given by $y(\mathbb{X}_{i})$, and QuAN is optimized through binary cross
entropy loss $\mathcal{L}=-\sum_{i}\hat{y_{i}}\log y(\mathbb{X}_{i})$ between
true label $\hat{y_{i}}$ and output.

###  A.2 Convolution

Figure 1: Schematic example of convolution layer applied to $N_{q}=36$ bit-
string with $N_{r}=6$ rows and $N_{c}=6$ columns. Convolution operation
denoted as ‘$*$’ involves summing over element-wise multiplication of part of
bit-string and convolution filter. $n_{c}$ represents the number of $2\times
2$ convolution filter, and $\mathbb{x}^{\alpha}$ is the output after
flattening convolution output into vector. We use a stride of $1$ and no
padding.

The original input data are sets of binary-valued arrays. We first pass the
input sets through a convolution layer (see SFig. 1). The convolution step has
two purposes: First, the original binary-valued arrays are mapped to vectors
with continuous entries, which have better algebraic properties. This is
illustrated in SFig. 1. Second, the convolution enables the model to capture
possible local spatial features. In the convolution layer, we apply
convolution filter $\\{F^{c}\\}_{c=1}^{n_{c}}$ of kernel size
$\texttt{kernel}=2$ and stride 1, on each two-dimensional bit-string array of
$N_{r}$ rows and $N_{c}$ columns. The channel number $n_{c}$ is a
hyperparameter of the ML model that controls the numbers of $2\times 2$
filters. After the convolution layer, BatchNorm2d over set elements follows.
The resulting output is then flattened into a 1D vector. The dimension of 1D
vector $d_{x}$ depends on the convolution layer hyperparameters;
$d_{x}=n_{c}(N_{r}-\texttt{kernel}+1)(N_{c}-\texttt{kernel}+1)$. We use a
stride of 1 and no padding in this paper. The output of the convolution layer
$\mathbb{x}^{\alpha}_{\mu}\equiv\mathbb{x}_{i,\mu}^{\alpha}$ is a matrix of
size $(N,d_{x})$, where $\alpha=1,\cdots,N$ and $\mu=1,\cdots,d_{x}$.

###  A.3 Mini-set Self-attention block (MSSAB) Encoder

The encoder aims to transform its input into a more informative
representation, which will be further used in the following decoder block to
accomplish the desired task — in our case, the binary classification task.
When designing the encoder architecture, we consider two key insights into the
data of interest: First, we expect the high-order moments of the bit-string
distribution to contain important information. Second, we expect the ordering
of the snapshots to be irrelevant since each measurement is done
independently.

The first insight, the desire to capture high-order moments, motivates us to
utilize the attention mechanism. The attention mechanism first introduced in
Ref. [1] drives the success of transformers as the core of large language
models. The attention mechanism lets the models learn correlations between
words (or tokens) in sentences. However, direct usage of the vanilla attention
mechanism brings us very limited power since the model will view the input
bit-strings as a sequence and will try to learn from their ordering.

To overcome this limitation, we need to consider the second insight, the
permutation invariance of the input bit-strings. Building on the vanilla
attention mechanism [1], Ref. [3] introduces a permutation invariant version
of transformer, named set transformer, since it treats the input as a set,
instead of a sequence. We follow the convention of the set transformer [3] and
call the self-attention module respecting permutation invariance the self-
attention block (SAB). Although, in theory, SAB can capture correlations while
respecting permutation invariance, a long sequence of stacked SABs is needed
for the model to access high-order moments, which soon becomes impractical.

Figure 2: Details of various blocks in the QuAN architecture. (a) The inner
structure of parallel self-attention block (SAB). (b) Multihead-attention
block (MAB) inside recurrent and reducing attention block (RecAB, RedAB). (c)
The inner structure of recurrent attention block (RecAB). (d) The inner
structure of reducing attention block (RedAB). $\sigma:S\to S$ is the random
permutation function that permutes mini-set index
$S=\\{0,1,\cdots,N_{s}-1\\}.$

We introduce MSSAB as a parameter-efficient and practical version of SAB when
accessing high-order inter-snapshot correlation between bit-strings. While the
SAB [3] considers all-to-all second-order inter-snapshot correlations between
input set elements, our MSSAB samples higher-order inter-snapshot correlations
and thus greatly reduces the computational cost. A single-layer of MSSAB
consists of the following three parts: parallel self-attention block (SAB),
recurrent attention block (RecAB), and reducing attention block (RedAB) (see
main text Fig. 1(e) encoder and SFig. 2). We discuss each part in more detail
in Sections A.3.1, A.3.2, A.3.3. Finally, we compare the computational
complexity of MSSAB with SAB in Section A.3.4.

####  A.3.1 Parallel SAB

First, the input set is shuffled and partitioned into $N_{s}$ mini-sets, where
$N_{s}$ is the number of mini-sets and is the core hyperparameter of MSSAB.
The parallel SAB transforms each mini-set independently. Note that the input
set to the encoder is the output of the previous convolution layer.

The parallel SAB block is shown as the purple block in main text Fig. 1(e),
and the inner structure is plotted in SFig. 2(a). We first perform a
preprocessing strategy called ‘mini-set partitioning’, which is a partition of
the input set (set size $N$) into $N_{s}$ subsets of size $N/N_{s}$, called
mini-sets. These mini-sets allow for parallel processing. The main objective
of parallel SAB is to capture pairwise second-order correlation only within
each mini-set.

The essential unit in parallel SAB is the self-attention block (SAB) that
allows access to second-order moments between bit-strings, which we will
discuss at the end of this paragraph. In the main text, we introduce the
simplified version of the self-attention score for a set of snapshots
$\mathbb{X}_{i}$ (see main text Eq. (1)). For actual implementation, we use
the output of convolution layer $\mathbb{x}_{i}=F^{c}(\mathbb{X}_{i})$ instead
of bare snapshots,

$\langle
Q\mathbb{x}_{i}|K\mathbb{x}_{i}\rangle=(Q\mathbb{x}_{i})(K\mathbb{x}_{i})^{T}$
(1)

up to normalization. We will omit index $i$ for convenience. The input to the
encoder is denoted as $\mathbb{x}_{\mu}^{\alpha}$ where $\alpha$ indexes the
set element and $\mu$ represents the dimension of the feature space after the
previous convolution layer. We will use the Greek letter $\alpha,\beta$ to
represent a set index running from $1$ to $N$ and $\mu,\nu,\rho,\lambda,\eta$
for the feature space index running from $1$ to $d_{h}$ (or $d_{x}$). SAB
transforms the input set to a set of hidden state vectors:

$\mathbb{h}_{\mu}^{\alpha}=\sum_{\nu}^{d_{x}}Q_{\mu\nu}\mathbb{x}_{\nu}^{\alpha}+\sum_{\beta=1}^{N}\texttt{Softmax}\left[\sum_{\rho}^{d_{h}}\sum_{\lambda\eta}^{d_{x}}\frac{1}{\sqrt{d_{h}}}\left(Q_{\rho\lambda}\mathbb{x}_{\lambda}^{\alpha}K_{\rho\eta}\mathbb{x}_{\eta}^{\beta}\right)\right]\sum_{\nu}^{d_{x}}V_{\mu\nu}\mathbb{x}^{\beta}_{\nu},$
(2)

with the hidden dimension size $d_{h}$. $Q$, $K$ and $V$ are query, key, and
value matrices of dimensions $(d_{h},d_{x})$. To gain insight into how the
model learns important relevant features, we rewrite the above expression in
the following form:

$\displaystyle A^{\alpha\beta}$
$\displaystyle=\texttt{Softmax}\left[\sum_{\rho\lambda\eta}\frac{1}{\sqrt{d_{h}}}\left(Q_{\rho\lambda}\mathbb{x}_{\lambda}^{\alpha}K_{\rho\eta}\mathbb{x}_{\eta}^{\beta}\right)\right]\sim\texttt{Softmax}\left[\mathbb{x}^{\alpha}\cdot\mathbb{x}^{\beta}\right],$
(3) $\displaystyle\mathbb{h}_{\mu}^{\alpha}$
$\displaystyle=\sum_{\beta}A^{\alpha\beta}\left(\sum_{\nu}Q_{\mu\nu}\mathbb{x}_{\nu}^{\alpha}+V_{\mu\nu}\mathbb{x}_{j}^{\beta}\right),$
(4)

where $\sum_{\beta}A^{\alpha\beta}=1$. The self-attention score matrix
$A^{\alpha\beta}$ is of dimensions $(N,N)$. Calculating the self-attention
score involves two set elements ($\mathbb{x}^{\alpha}$ and
$\mathbb{x}^{\beta}$), which can capture all-to-all second-order moments in
$\mathbb{x}$. This is followed by layer normalization on spatial dimension and
a linear layer:

$\displaystyle\mathbb{h}_{\mu}^{\prime\alpha}$
$\displaystyle=\texttt{LayerNorm}(\mathbb{h}_{\mu}^{\alpha}),$ (5)
$\displaystyle\mathbb{y}_{\mu}^{\alpha}$
$\displaystyle=\texttt{Sigmoid}\left(\texttt{LayerNorm}(\mathbb{h}^{\prime}_{\mu}+\texttt{FF}_{\mu\nu}(\mathbb{h}^{\prime}_{\nu}))\right),$

where FF is a feed-forward function for residual connection that acts on each
set element equally, which in our implementation is to multiply by a matrix
$O$ of dimensions $(d_{h},d_{h})$ followed by activation function (either
Sigmoid or ReLu)), i.e. $\texttt{FF}(\mathbb{x})=\texttt{ReLu}(O\mathbb{x})$
or $\texttt{Sigmoid}(O\mathbb{x})$ for each set element $\mathbb{x}$.
$Q,K,V,O$ are learned weight matrices in SAB. The output of SAB,
$\mathbb{y}^{\alpha}_{i}$, is a matrix of dimensions $(N,d_{h})$.

There are $N_{s}$ parallel SABs acting on $N_{s}$ mini-set, where parameters
are shared across all the SAB blocks. We denote the output of parallel SAB as
$\\{\mathbb{y}_{(0)},\cdots,\mathbb{y}_{(N_{s}-1)}\\}$, where the subscript
with parentheses is the mini-set label. Each output $\mathbb{y}_{(m)}$ is a
matrix of size $(N/N_{s},d_{h})$. Parallel SAB can access second-order moments
of a bit-string distribution within each mini-set, in contrast to ordinary SAB
accessing all-to-all second-order moments.

####  A.3.2 Recurrent AB (RecAB)

The RecAB takes the outputs of parallel SAB, and attends them recurrently with
multiple randomized orderings. A schematic of RecAB is shown as the blue block
in Fig. 1(e), and the inner structure is plotted in SFig. 2(c). RecAB was
devised to capture correlations between different mini-sets that were not
captured through parallel SAB; by attending mini-sets
$\\{\mathbb{y}_{(0)},\cdots,\mathbb{y}_{(N_{s}-1)}\\}$ recurrently, RecAB can
capture up to $2N_{s}$-th order correlation in $\mathbb{x}$.

We utilize multihead-attention block (MAB) instead of SAB since we compute the
attention score between two different mini-sets. Each (parallel SAB output)
mini-set $\mathbb{y}_{(m)}$ goes into $N_{s}-1$ MABs. At each time $t$
attended by the next mini-set $\mathbb{y}_{((m+t+1)\texttt{mod}N_{s})}$, where
MAB is given by

$\mathbb{h}_{(t+1),\mu}^{\alpha}=\sum_{\nu}^{d_{h}}Q^{\prime}_{\mu\nu}\mathbb{y}_{((m+t+1)\texttt{mod}N_{s}),\nu}^{\alpha}+\sum_{\beta=1}^{N/N_{s}}\texttt{Softmax}\left[\sum_{\rho\lambda\eta}^{d_{h}}\frac{1}{\sqrt{d_{h}}}\left(Q^{\prime}_{\rho\lambda}\mathbb{y}_{((m+t+1)\texttt{mod}N_{s}),\lambda}^{\alpha}K^{\prime}_{\rho\eta}\mathbb{h}_{(t),\eta}^{\beta}\right)\right]\sum_{\nu}^{d_{h}}V^{\prime}_{\mu\nu}\mathbb{h}_{(t),\nu}^{\beta},$
(6)

where $\mathbb{h}_{(0)}\equiv\mathbb{y}_{(m)}$,
$\mathbb{h}_{(N_{s}-1)}\equiv\mathbb{y}^{\prime}_{(m)}$. Attention score
between two different mini-sets involves two set elements
$\mathbb{y}^{\alpha}_{(m)}$ and $\mathbb{y}^{\beta}_{(m^{\prime})}$ each from
$m$ and $m^{\prime}$-th mini-set, hence capturing the second-order correlation
between two mini-sets. Each MAB operation is independent and identical,
followed by the same layers as in Eq. (5). We denote the output of RecAB as
$\\{\mathbb{y}_{(0)}^{\prime},\cdots,\mathbb{y}_{(N_{s}-1)}^{\prime}\\}$. Each
output $\mathbb{y}_{(m)}^{\prime}$ is a matrix of size $(N/N_{s},d_{h})$.
While each mini-set passes through $N_{s}-1$ numbers of MAB recurrently, it
involves $N_{s}$ set elements each from different mini-set
($\mathbb{y}^{\alpha_{0}}_{(m)},\mathbb{y}^{\alpha_{1}}_{(m+1)},\cdots,\mathbb{y}^{\alpha_{N_{s}-1}}_{((m+N_{s}-1)\texttt{mod}N_{s})}$),
which can capture $N_{s}$-th order moment in $\mathbb{y}$. In other words,
RecAB can access $2N_{s}$-th order moments in $\mathbb{x}$, considering
$\mathbb{y}$ contains second-order moment information of $\mathbb{x}$.

####  A.3.3 Reducing AB (RedAB)

Finally, the RedAB attends mini-sets in a randomized sequence and shrinks
$N_{s}$ mini-sets into one.

The RedAB is shown as the blue dashed block in Fig. 1(e), and the inner
structure is plotted in SFig. 2(d). RedAB is designed to reduce all (RecAB
output) mini-sets into a single mini-set while preserving the mini-set
permutation invariance. We attend mini-sets
$\\{\mathbb{y}_{(0)}^{\prime},\cdots,\mathbb{y}_{(N_{s}-1)}^{\prime}\\}$ in a
randomized sequence using $N_{s}-1$ MABs (each MAB operation is independent
and identical, followed by the same layers as in Eq. (5); we use the same MAB
from RecAB.) Similar to RecAB, RedAB with $N_{s}-1$ numbers of MAB involves
$N_{s}$ set elements from different mini-set
($\mathbb{y}^{\prime\beta_{0}}_{(\sigma(0))},\mathbb{y}^{\prime\beta_{1}}_{(\sigma(1))},\cdots,\mathbb{y}^{\prime\beta_{N_{s}-1}}_{(\sigma(N_{s}-1))}$
where $\sigma$ is a randomized permutation function), which captures
$N_{s}$-th order moment in $\mathbb{y}^{\prime}$. In other words, RedAB can
access $2N_{s}^{2}$-th order moments in $\mathbb{x}$. The final output of
RedAB is $\mathbb{z}\equiv\mathbb{z}^{\alpha}_{\mu}$, a matrix of dimensions
$(N/N_{s},d_{h})$ with $\alpha=1,\cdots,N/N_{s}$ and $\mu=1,\cdots,d_{h}$.
Therefore, one MSSAB layer reduces the input set size $N$ to $N/N_{s}$.

####  A.3.4 Computational complexity of MSSAB

MSSAB is more parameter-efficient when we need to target high-order moments.
In this subsection, we discuss the computational complexity of MSSAB and
compare it to the SAB [3].

First, we walk through how the MSSAB collects increasing order of moments
through the parallel-SAB, RecSAB, and RedSAB sequence and finally reaches up
to $(2N_{s}^{2})$-th moments in its output.

Parallel-SAB -  In parallel-SAB block, each mini-set passes through one layer
of SAB (see SFig. 2(a)). The transformation performed in SAB is shown in Eq.
(2), where two set elements are involved in the calculation of attention score
(see Eq. (3)), so the outputs $\\{\mathbb{y}_{(m)}\\}_{m=0}^{N_{s}-1}$ samples
up to $2$-nd order moments of the input set $\mathbb{x}$. For simplicity, we
define $\texttt{Order}(\cdot)$ to be the order of moments $(\cdot)$ can
access. Thus, $\mathtt{Order}(\mathbb{y}_{(m)})=2$ for
$\mathbb{y}_{(m)}\in\\{\mathbb{y}_{(m)}\\}_{m=0}^{N_{s}-1}$.

RecAB -  The RecAB attends each of $\\{\mathbb{y}_{(m)}\\}_{m=0}^{N_{s}-1}$
with others, resulting in $\\{\mathbb{y}_{(m)}^{\prime}\\}_{m=0}^{N_{s}-1}$.
It is more convenient to calculate the order of moments if we consider the
recurrent representation, as shown in SFig. 2(b). Recursively, we have

$\displaystyle\mathtt{Order}(\mathbb{h}_{(t+1)})$
$\displaystyle=\mathtt{Order}(\mathbb{h}_{(t)})+\mathtt{Order}(\mathbb{y}_{(t+1)\text{mod}N_{s}})$
(7) $\displaystyle=\mathtt{Order}(\mathbb{h}_{(t)})+2,$ (8)
$\displaystyle\mathtt{Order}(\mathbb{h}_{(0)})$
$\displaystyle=\mathtt{Order}(\mathbb{y}_{(m)})=2.$ (9)

Via recursion unrolling, we have
$\mathtt{Order}(\mathbb{h}_{N_{s}-1})=2N_{s}$. As shown in the unrolled
flowchart of RecAB (SFig. 2(c)), each row represents a different ordering of
$\\{\mathbb{y}_{(m)}\\}_{m=0}^{N_{s}-1}$ in which the attention mechanism
attends. That is, we initialize $\mathbb{h}_{(0)}$ with $N_{s}$ different
mini-sets, each gives an output among
$\\{\mathbb{y}_{(m)}^{\prime}\\}_{m=0}^{N_{s}-1}$. Thus,

$\mathtt{Order}(\mathbb{y}_{(m)}^{\prime})=\mathtt{Order}(\mathbb{h}_{N_{s}-1})=2N_{s}.$
(10)

RedAB -  RedAB is similar to the RecAB since both use a recurrent module.
However, RedAB only samples one order in which the outputs of RecAB
$\\{\mathbb{y}_{(m)}^{\prime}\\}_{m=0}^{N_{s}-1}$ get attended. Again, using
the recurrent representation, we now have

$\displaystyle\mathtt{Order}(\mathbb{h}^{\prime}_{(t+1)})$
$\displaystyle=\mathtt{Order}(\mathbb{h}^{\prime}_{(t)})+\mathtt{Order}(\mathbb{y}^{\prime}_{(\sigma(t+1)})$
(11) $\displaystyle=\mathtt{Order}(\mathbb{h}^{\prime}_{(t)})+2N_{s},$ (12)
$\displaystyle\mathtt{Order}(\mathbb{h}^{\prime}_{(0)})$
$\displaystyle=\mathtt{Order}(\mathbb{y}_{\sigma(0)})=2N_{s}.$ (13)

Unrolling the recursive relation, we obtain

$\mathtt{Order}(\mathbb{z})=\mathtt{Order}(\mathbb{h}^{\prime}_{N_{s}-1})=2N_{s}^{2}.$
(14)

We have shown that a single layer of MSSAB can access moments up to
$2N_{s}^{2}$-th order. By stacking $L$ layers of MSSAB, we can reach
$(2N_{s}^{2})^{L}$-th moments.

Second, we consider the runtime complexity of the MSSAB versus the SAB. In a
SAB, attention scores are computed between each pair of set elements,
accessing up to second-order moments. Consequently, the computational
complexity scales quadratically with the input set size $N$, resulting in an
$O(N^{2})$ complexity. Our MSSAB remains a similar quadratic scaling of
complexity $O(N^{2})$. Therefore, accessing higher moments with much fewer
layers makes MSSAB more efficient than SAB.

Figure 3: Number of layers $L$ required for SAB and MSSAB to reach moments of
order $\Theta$. For simplicity, we omit the ceiling ($\lceil\cdot\rceil$)
operation in the plot. For SAB, the green curve plots $L=\log_{2}\Theta$. For
MSSAB, we plot $L=\log_{2N_{s}^{2}}\Theta$ for $N_{s}=5$ (red curve). The red
star marks 1 layer of MSSAB with $N_{s}$, which is the setup used in our RQC
training in SM section D.

To reach moments of order $\Theta$, one need
$L=\lceil\log_{2N_{s}^{2}}\Theta\rceil$ layers of MSSAB. In contrast, reaching
the same order needs $\lceil\log_{2}\Theta\rceil$. We show in SFig. 3 a
scaling of the number of layers required for SAB or MSSAB to reach moments of
desired order $\Theta$. Clearly, to reach larger $\Theta$, the number of
layers needed for SAB grows much faster than MSSAB. With $L=1$ layer of MSSAB
with $N_{s}=5$, we can access $\Theta=50$, which is marked by a red star in
SFig. 3. In contrast, we need $L=6$ SAB layers to reach the same order
$\Theta=50$ with SAB.

###  A.4 Pooling Attention Block (PAB) Decoder

While the MSSAB encoder handles various orders of moments between bit-strings,
the main operation of the decoder is to pool out useful information from the
encoder output while respecting the permutation invariance of the snapshots.
Unlike ordinary pooling operations (such as averaging/summing over snapshots),
pooling attention block (PAB) enables importance-sampling by assigning weight
through pooling attention score.

The decoder of QuAN consists of one layer of pooling attention block (PAB) and
a final single-layer perceptron to output a single scalar value, which is the
label prediction confidence. The output of the encoder
$\mathbb{z}_{\mu}^{\alpha}$ is fed into the decoder as its input. Note that so
far, the set dimension, indexed by $\alpha$, has the dimensionality as the
original input $\mathbb{x}$ divided by $N_{s}$. In the decoder, the essential
operation is done by the PAB, which weighs different set elements differently
so that elements with important features contribute more to the final result.
Writing explicitly, the PAB first transforms the encoder output as

$\displaystyle\mathbb{p}_{\mu}$
$\displaystyle=S_{\mu}+\sum_{\beta=1}^{N/N_{s}}\texttt{Softmax}\left[\sum_{\rho\lambda}^{d_{h}}\frac{1}{\sqrt{d_{h}}}\left(S_{\rho}K^{\prime\prime}_{\rho\lambda}\mathbb{z}_{\lambda}^{\beta}\right)\right]\sum_{\nu}^{d_{h}}V^{\prime\prime}_{\mu\nu}\mathbb{z}_{\nu}^{\beta}$
(15)
$\displaystyle=\sum_{\beta}s^{\prime\beta}\left(S_{\mu}+\sum_{\nu}V^{\prime\prime}_{\mu\nu}\mathbb{z}_{\nu}^{\beta}\right),$

where a seed vector $S$ is used as query vector of size $(1,d_{h})$ for a
weighted average $\mathbb{z}^{\beta}$ over the set dimension, and
$K^{\prime\prime}$, $V^{\prime\prime}$ are key and value matrices of
dimensions $(d_{h},d_{h})$. This operation plays an important role in making
the output $y(\mathbb{X})$ set permutation invariant, which means the output
is the same even when shuffling the set element. Moreover, the pooling
attention score $s^{\prime\beta}$ can be considered as a weight of
$\mathbb{z}^{\beta}$; therefore, $\mathbb{p}$ can be considered as a weighted
sum of encoder output. (See SM section E.5 for a detailed discussion on
pooling attention score in identifying topologically ordered quantum states.)
Here, the pooling attention score $s^{\prime\beta}$ is a vector of size $N$:

$s^{\prime\beta}=\texttt{Softmax}\left[\sum_{\rho\lambda}\frac{1}{\sqrt{d_{h}}}\left({S}_{\rho}K^{\prime\prime}_{\rho\lambda}\mathbb{z}_{\lambda}^{\beta}\right)\right].$
(16)

Similar to SAB in the encoder, we perform a layer normalization on spatial
dimension and a linear layer:

$\displaystyle\mathbb{p}_{\mu}^{\prime}$
$\displaystyle=\texttt{LayerNorm}(\mathbb{p}_{\mu})$ (17) $\displaystyle
y(\mathbb{X})$
$\displaystyle=\texttt{Sigmoid}\left(\sum_{\mu}W_{\mu}\;\texttt{LayerNorm}\left[\mathbb{p}_{\mu}^{\prime}+\texttt{rFF}_{\mu\nu}(\mathbb{p}^{\prime}_{\nu})\right]+b\right)$
(18)

where $W$ is a matrix of dimensions $(1,d_{h})$ that converts the vector
$\mathbb{p}_{\mu}^{\prime}$ into scalar output $y$, which is the confidence of
predicting given set $\mathbb{X}$ into one of the classes (e.g. volume-law
phase, deep circuit outcome or topological phase for the cases in the main
text). $S,K^{\prime\prime},V^{\prime\prime},W,b$ are learnable parameters in
the decoder.

##  Appendix B Comparison of QuAN and other ML architectures

In this section, we provide details of all the ML architectures that are
compared with QuAN in Fig. 1(e). In particular, we highlight the main
differences between these architectures and the QuAN architecture we propose
in this work. The models maintain approximately the same total number of
trainable parameters to make a controlled comparison between different
architectures.

##### Multi-layer perceptron (MLP) -

The multi-layer perceptron (MLP) model contains 4 layers of perceptrons in the
encoder, followed by the decoder block composed of one linear layer and
another perception (see SFig. 4). The first linear layer in the decoder
transforms each set element independently, while the second linear layer
combines information from different set elements. Note that the MLP
architecture does not use the attention mechanism and also does not respect
set element permutation invariance.

Figure 4: Schematic MLP architecture. The inner structure of the perceptron
block (marked as the orange box) is shown in the right panel. The dash-lined
blocks indicate inter-set-element operations. Blocks with solid lines act on
each set element independently.

##### Convolutional neural network (CNN) -

The convolutional neural network (CNN) architecture (see SFig. 5) treats the
input data as a 3-dimensional object, hence it does not respect the
permutation invariance of set elements either. The input 3D array goes over 3
layers of normalized 3D convolution that extracts relevant features into
$n_{c}$ channels, and the subsequent linear layer and sigmoid function outputs
a binary classification probability.

Figure 5: Schematic CNN architecture. The dash-lined blocks indicate inter-
set-element operations. Blocks with solid lines act on each set element
independently.

##### Transformer (Transf.) -

The transformer attends to spatial correlation within one measurement outcome.
In contrast to all the other models, the transformer takes one bit-string as
an input at a time, considering each bit (0 or 1) as a token. The bit-string
goes through a 2D positional encoding followed by ordinary self-attention
blocks, a linear layer, and an activation function, as shown in SFig. 6.
Compared to ordinary transformers, our Transf. lacks a decoder part, i.e.,
this model is analogous to an ordinary transformer encoder.

Figure 6: Schematic transformer architecture.

##### Set multi-layer perceptron (SMLP) -

The set multi-layer perceptron (SMLP) is similar to the MLP architecture
defined above. The main difference between MLP and SMLP lies in the decoder,
shown in SFig. 7. For SMLP, we modify the decoder such that the model respects
the set element permutation invariance. SMLP only contains linear layers,
summation pooling, and final non-linear activation functions. The attention
mechanism is absent in this model.

Figure 7: Schematic SMLP architecture. The dash-lined blocks indicate inter-
set-element operations. Blocks with solid lines act on each set element
independently.

##### Pooling attention block (PAB) -

In PAB architecture shown in SFig. 8., we replace the encoder of QuAN with two
MLP layers. The decoder still contains the PAB block, which utilizes the
attention mechanism.

Figure 8: Schematic PAB architecture. The dash-lined blocks indicate inter-
set-element operations. Blocks with solid lines act on each set element
independently.

##  Appendix C Driven Hard-core Bose-Hubbard model

###  C.1 Data acquisition

In this section, we provide further details of learning the entanglement
transition in the driven hard-core Bose-Hubbard model.

We use the data acquired from a $4\times 4$ array of superconducting transmon
qubits as described in Ref. [4]. In this system, the on-site interaction,
determined by the anharmonicity of the transmon qubits, is much stronger than
the exchange interaction. Therefore, the system can be described by a Hard-
core Bose-Hubbard (HCBH) Hamiltonian:

$\hat{H}_{\mathrm{HCBH}}/\hbar=\sum_{\langle
j,k\rangle}J_{jk}\hat{\sigma}_{j}^{+}\hat{\sigma}_{k}^{-}-\sum_{j}\frac{\epsilon_{j}}{2}\hat{\sigma}^{z}_{j}\,,$
(19)

where $\hat{\sigma}_{j}^{+}$ ($\hat{\sigma}_{j}^{-}$) is the raising
(lowering) operator for a two-level system at site $j$, and
$\hat{\sigma}^{z}_{j}$ is the Pauli-$Z$ operator. The first term describes the
particle exchange interaction between neighboring lattice sites with strength
$J_{jk}$, with an average strength of
$J/2\pi=$5.9(4)\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$. The second term
represents the site energies, which are tuned by the transmon transition
frequencies with an accuracy of $300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$
($\approx 5\times 10^{-2}J$) in the device [5]. This system features site-
resolved, multiplexed single-shot dispersive qubit readout, which, combined
with single-qubit gates, enables simultaneous tomographic measurements of the
qubit states. In order to prepare superposition states across the energy
spectrum of the lattice, the interacting qubits are simultaneously driven via
a common control line. The Hamiltonian of the driven lattice is

$\hat{H}/\hbar=\sum_{\langle
j,k\rangle}J_{jk}\hat{\sigma}_{j}^{+}\hat{\sigma}_{k}^{-}+\frac{\delta}{2}\sum_{j}\hat{\sigma}^{z}_{j}+\Omega\sum_{j}(\alpha_{j}\hat{\sigma}_{j}^{-}+\mathrm{h.c.})\,,$
(20)

where $\delta$ is the detuning between the drive and the qubit frequencies
(all sites have the same energy). The drive strength $\Omega$ can be tuned by
varying the amplitude of the applied drive pulse. The common drive couples to
each qubit with a complex coefficient $\alpha_{k}$ (see Ref. [4] for details).

By changing the drive detuning $\delta$, the distribution of the superposition
states across the HCBH energy spectrum can be controlled: with detuning
$\delta=0$, the superposition state will be concentrated near the center of
the energy band, whereas as the magnitude of $\delta$ increases, the
superposition state approaches the edge of the energy band. Therefore, the
drive detuning is an effective tuning knob to control the distribution of the
state across the spectrum to study the entanglement scaling behavior in many-
body systems.

###  C.2 Data preprocessing

The states $\rho=\rho(\delta,t)$ accessible in the emulation is parameterized
by the driving detuning parameter $\delta$ and the driving time $t$. For each
given state $\rho(\delta,t)$, we acquire a corresponding probability
distribution $p(b)=\text{Tr}(\rho|b\rangle\langle b|)$ by experimentally
measuring the state $10^{4}$ times in the $Z$-basis. We then sample snapshots
from each state’s probability $p(b)$ to generate training and testing
datasets. Given that the increase in detuning strength $\delta$ correlates
directly with the total particle number, we selectively remove snapshots
featuring a total particle number other than $n=8$. Such filtering prevents
the ML model from learning the trivial feature – the particle density. The
resulting number of snapshots we keep for each state is $M_{s}=4096$, each
having $n=8$ total particle number. We emulate $9$ different values of
$\delta/J$, namely $\delta/J\in\\{-2,-1.5,\cdots,0,\cdots,1.5,2\\}$. For each
fixed value of $\delta/J$, we take $17$ grid points along $t\in\\{1.4\times
10^{-7},1.5\times 10^{-7},\cdots,2.9\times 10^{-7},3\times 10^{-7}\\}$
seconds, all of which are in the steady state region. Since we are interested
in steady states at different $\delta$, we batch along $t$. Hence, for each
$\delta$ we have $M=4096\times 17=69632$ snapshots. Our device is a $4\times
4$ array of superconducting transmon qubits, so each snapshot also has a
$4\times 4$ rectangular geometry with $N_{q}=16$ binary values.

We train our model on snapshots taken from the volume-law states at
$\delta/J=0$ and from the area-law states at $\delta/J=\pm 2$. For each
$\delta/J$, we randomly choose 12 out of 17 different values of $t$ for
training and keep the remaining 5 for validation. We then randomly partition
the $M_{s}=4096$ snapshots from each state $\rho(\delta,t)$ into sets of size
$N$. Half of the sets from $\delta=\pm 2$ are combined and used in training so
that the total training data size remains the same for the two classes. To
train QuAN as a binary classifier between the two phases, we perform
supervised learning by labeling sets from $\rho(\delta/J=0,t)$ as $\hat{y}=0$,
and $\rho(\delta/J=\pm 2,t)$ as $\hat{y}=1$. The details of the training
procedure are discussed in the next section.

After training the model, we test the model on testing dataset taken from the
intermediate phase region $\delta/J\in\\{-1.5,\cdots,0,\cdots,1.5\\}$. The
testing sets are generated by the same procedure as the training sets but at
different detuning strength $\delta/J$. We obtain average confidence $\bar{y}$
over testing sets in predicting the volume-law phase. For $\delta/J=0$ and
$\delta/J=\pm 2$, we present the average confidence for the reserved
validation sets.

Model hyperparameter
---
Model | SMLP, PAB, QuAN2, QuAN4
Set size ($N$) | 1$\sim$256
Number of mini-sets ($N_{s}$) | 1
Number of MSSAB layer ($L$) | 0$\sim$2
Number of $2\times 2$ Conv. channel ($n_{c}$) | 7$\sim$8
Attention block: Hidden spatial dimension ($d_{h}$) | $16$
Attention block: Number of heads ($n_{h}$) | $4$
Attention block: activation function for residual connection | Sigmoid
Training hyperparameter
Optimizer |  `Adam`$(\beta_{1}=0.9,\;\beta_{2}=0.999,\;\epsilon=1\times 10^{-8})$
L2 coefficient | $5\times 10^{-5}$
Learning rate | $1\times 10^{-4}$
LR schedule |  `StepLR`$(\text{stepsize}=200,\gamma=0.65)$
Epoch | 500
Dataset shuffling period | $10$
Batchsize | $80000/N$
Initialization | Default
GPU | A100 (80GB)
Table 1: Model setting and training hyperparameters for the driven hard-core Bose-Hubbard model. | SMLP | PAB | QuAN2 | QuAN4
---|---|---|---|---
| Conv(1, 8, 2, 1, BatchNorm) | Conv(1, 8, 2, 1, BatchNorm) | Conv(1, 8, 2, 1, BatchNorm) | Conv(1, 7, 2, 1, BatchNorm)
Encoder | SLP(72, 48, Sigmoid) | SLP(72, 48, Sigmoid) | MSSAB*(72, 16, 4, 1) | MSSAB*(63, 16, 4, 1)
| SLP(48, 16, Sigmoid) | SLP(48, 16, Sigmoid) |  | MSSAB*(16, 16, 4, 1)
| SLP(16, 48, Sigmoid) | PAB*(16, 16, 4) | PAB*(16, 16, 4) | PAB*(16, 16, 4)
Decoder | SLP(48, 1) | SLP(16, 1, Sigmoid) | SLP(16, 1, Sigmoid) | SLP(16, 1, Sigmoid)
| Sum*, Sigmoid |  |  |
Table 2:  Architecture parameters used in models sketched in main text Fig.
2(e-g). Data from $N_{q}=16$ with $4\times 4$ geometry was used. Each layer’s
arguments are as follows: Convolutional layer as Conv($n_{c,\text{in}}$,
$n_{c,\text{out}}$, kernel, stride, normalization). Single-layer perception as
SLP($d_{\text{in}}$, $d_{\text{out}}$, activation). MSSAB($d_{h,\text{in}}$,
$d_{h,\text{out}}$, $n_{h}$, $N_{s}$). PAB($d_{h,\text{in}}$,
$d_{h,\text{out}}$, $n_{h}$). The asterisk (*) denotes the module operates on
a set dimension.

###  C.3 Training and testing procedure

We use PyTorch to train and test the model to discriminate volume-law and
area-law entanglement scalings. Using the Adam optimization algorithm, we
minimize the binary cross entropy loss function (`torch.nn.BCEloss`) of the
true output and machine-predicted output of the given input set. To prevent
the model from overfitting, we employ the “dataset shuffling period”: since
the model input is set structured with multiple measurement outcomes, we
shuffle or regenerate the input set-structured dataset every 10 epochs to
ensure the model can explore various combinations of measurement outcomes.
This scheme is used throughout the paper. Training hyperparameters are listed
in Table 1, and parameters of architectures are listed in Table 2. We train
different models and different set sizes (ranging from $N=2^{0}=1$ to
$N=2^{8}=256$) independently while keeping the remaining model parameters
unchanged. For each architecture, we perform $10$ independent training to
ensure the stability of training. We store the model with the highest accuracy
on validation data. When calculating accuracies, we set the classification
criteria (threshold) to be $y=0.5$; that is, $\mathbb{X}_{i}$ is classified as
volume-law if the machine output is $y(\mathbb{X}_{i})>0.5$ and is classified
as area-law otherwise. Comparing the machine-predicted labels for validation
sets with the expected labels yields the validation accuracy.

We don’t know the labels in the intermediate region, where testing sets are
generated. However, we can see the phase transition from the machine
confidence $y(\mathbb{X}_{i})$. We get averaged confidence by averaging the
machine output of confidence $y(\mathbb{X}_{i})$ over testing sets. We test
$10$ independent models and obtain the mean and the standard error of the mean
of the machine confidences.

###  C.4 Machine learning details

Figure 9: Validation classification accuracy learning curves as a function of
training epoch. (a) Validation classification accuracy curves at different set
sizes $N$ using QuAN4. (b) Validation classification accuracy curves with
different architectures (SMLP, PAB, QuAN2 and QuAN4) using $N=64$. The solid
line shows the mean of 10 independently trained models, and shaded regions
show the minimum and maximum accuracy at each epoch.

After training, we verify the model’s performance and stability using
validation accuracy as a function of epoch (see SFig. 9). Validation accuracy
is determined by the percentage of correctly classified validation data points
that are in either volume-law ($\delta/J=0$) or area-law ($\delta/J=\pm 2$)
state (see section C.2). For instance, an accuracy $90\%$ indicates that the
model correctly classifies $90\%$ of the total data. If the accuracy hovers
around $50\%$, it suggests that the model fails to train, as random guessing
would yield the same accuracy of $50\%$ for any binary classification task.

From SFig. 9(a), we observe the differences in accuracy curves with different
set sizes $N$ we use to train QuAN4. When the set size is small (e.g. $N<8$),
test accuracy saturates at a much lower value compared to the larger set size
($N\geq 64$). The saturated accuracy increases as we increase the set size up
to $N=64$, suggesting that the model can better classify the data if using
larger set sizes. From the trend as we increase the set size, it is clear that
having a set structure for input is essential for QuAN to learn the
entanglement transition. However, we observe accuracy is upper-bounded by
increasing the set size beyond $N>64$. For the two largest set sizes ($N=128$
and $256$), the accuracy does not go higher than the saturated accuracy for
$N=64$. Moreover, the variance of learning curves across $10$ independent
training gets larger with increasing set size $N$ due to a decrease in total
training points. Even the set structure plays an important role, using large
set size is limited by the total number of snapshots.

SFig. 9(b) shows the accuracy curves for different architectures. All 4
different models require input with a set structure, and we maintain the same
set size $N=64$ for fair comparison. Models with self-attention (QuAN2, QuAN4)
exhibit higher accuracy saturation than models without self-attention (PAB,
SMLP), which hovers around 60%. Especially, QuAN4 accuracy saturates faster
and higher compared to QuAN2, demonstrating having more SAB layer is
advantageous. Another thing to notice is that the range or variance of
accuracy across $10$ independent training runs (shaded region) varies among
architectures. Notably, QuAN4 demonstrates good stability with smaller
variance at later epoch$\approx 500$, compared to other architectures.

Figure 10: Average confidence $\bar{y}$ as a function of detuning $\delta/J$.
(a) Average confidence by QuAN4 with varying set sizes $N$. (b) Average
confidence with different architectures (SMLP, PAB, QuAN2 and QuAN4) using
$N=64$. The red stars indicate the training points $\delta/J=0,\pm 2$.

Following model validation, we now test the model performance in SFig. 10 by
examining the averaged confidence $\bar{y}$ averaged over the testing dataset
that includes intermediate region
$\delta/J\in\\{-1.5,-1.0,-0.5,0.5,1.0,1.5\\}$ (see also main text Fig. 2(e-g)
in the main text). If the averaged confidence $\bar{y}$ is near $0.5$, we
consider the model fails to distinguish volume-law and area-law features. In
SFig. 10(a), we present the average confidence for predicting volume-law,
using different set sizes $N$ for our model QuAN4. Notably, as we increase the
set size up to $N=64$, the transition prediction becomes sharper, emphasizing
the importance of utilizing a set structure. However, for larger set sizes
beyond $N>64$ (i.e. $N=128$ and $256$), the average confidence trend over
drive detuning $\delta/J$ remains the same with the $N=64$ training. Moreover,
the variance (error bar) of average confidence $\bar{y}$ across 10 independent
training increases with set size $N>64$, which aligns well with the accuracy
curve in SFig. 9(a). Despite the benefits of large set sizes, there is a
trade-off due to a fixed number of snapshots $M=69632$; by increasing the set
size, we sacrifice the total number of sets available. Specifically, we have a
total $M/N=1088$ sets for $N=64$, $M/N=544$ sets for $N=128$ and $M/N=272$
sets for $N=256$. Beyond a set size of $N=64$, we encounter a risk of model
overfitting to the training dataset if we don’t use a sufficient amount of
training data points. We conclude that a set size beyond $N=64$ does not have
an advantage in observing entanglement transition.

In SFig. 10(b), we present the entanglement transition witnessed with
different architectures while maintaining the same set size $N=64$. Models
with self-attention (QuAN2 and QuAN4) exhibit a sharper distinction between
area-law and volume-law phases compared to models without self-attention (PAB,
SMLP). Moreover, we observe QuAN4 achieves the smallest error bar (variance
across 10 independent training) except for SMLP, which already demonstrates
poor performance in observing entanglement transition. This again highlights
the advantage of QuAN4 that can access high moments between sampled bit-
strings.

##  Appendix D Random quantum circuit

###  D.1 Data acquisition

####  D.1.1 Quantum processor details and experimental procedure

The random quantum circuit experiment was run on a Google Sycamore processor
composed of 70 frequency-tunable transmon qubits with tunable couplers. The
quantum processor used has a similar design to Ref. [6] and was carried out on
the same processor used in previous works, where typical coherence times,
readout errors, and single and two-qubit gate errors on this particular chip
can be found [7, 8]. The two-qubit gates used for this experiment are iSWAP-
like gates with an iSWAP angle $\theta\approx 0.5\pi$ and conditional phase
angle $\phi\approx 0.1\pi$ [8]. We collected data on rectangular subarrays of
$N_{q}=20,25,30$ and $36$ qubits (see SFig. 11) with variable circuit depth
$d=4,6,8,10,12,14,16,18$ and $20$. For every $N_{q}$ and $d$, we collected
data on 50 different random circuit instances. Each instance contains a
different sequence of single-qubit gates randomly chosen from gate set
$\\{\sqrt{X^{\pm 1}},\sqrt{Y^{\pm 1}},\sqrt{W^{\pm 1}},\sqrt{V^{\pm 1}}\\}$,
with $W=(X+Y)/\sqrt{2}$ and $V=(X-Y)/\sqrt{2}$. For each of the $N_{c}=50$
random circuit instances we performed $M_{s}=500,000$ ($M_{s}=2,000,000$)
$Z$-basis measurements for $N_{q}=20,25,30$ ($N_{q}=36$). Thus, a total of
$50\times M$ bit strings were collected for each $(N_{q},d)$ pair.

Figure 11: Layouts of Google Sycamore processor with 70 qubits (dark grey
circles). The subarrays used for system sizes $N_{q}=20,25,30,$ and $36$ are
marked in colored boxes.

####  D.1.2 Simulation and linear cross-entropy benchmarking (XEB)

In main text Fig. 3(c,f), we show the linear cross entropy benchmark (XEB)
$\mathcal{F}_{\text{XEB}}(N_{q},d)$ as a function of circuit depth $d$ for
different system sizes $N_{q}$. Here, we present how we obtain the data.

Parallelized over 8 NVIDIA A100 GPUs, we can simulate up to $N_{q}=36$ qubits
exactly. We use Cirq [9] to simulate the same random quantum circuit instances
as used in the experiment. For each circuit instance, we evolve an all-zero
product $|0\rangle^{\otimes N_{q}}$ with the circuit to get a state vector
$|\psi_{s}(N_{q},d)\rangle$, where $d\in[4,6,\cdots,20]$ represents the depth
of the circuit, and $s\in[1,2,\cdots,N_{c}]$ represents $N_{c}=50$ different
circuit instances. To simulate measuring a state in $Z$-basis, we sample from
the distribution given by $|\psi_{s}(N_{q},d)|^{2}$. Similar to the sample
size we have in experiments, we draw $M_{s}=500,000$ samples for
$N_{q}=20,25,30$ and $M_{s}=2,000,000$ for $N_{q}=36$. The linear XEB is
defined as,

$\mathcal{F}_{\text{XEB}}(N_{q},d,s)=2^{N_{q}}\langle
p(B_{i})\rangle_{i}-1=2^{N_{q}}\sum_{b\in(0,1)^{\otimes N_{q}}}p(b)^{2}-1$
(21)

where $p(B_{i})=|\langle\psi_{s}(N_{q},d)|B_{i}\rangle|^{2}$. From a finite
set of samples, we can get an estimate of $\mathcal{F}_{\text{XEB}}$,

$\mathcal{F}_{\text{XEB}}(N_{q},d,s)\approx
2^{N_{q}}\left(\frac{1}{M}\sum_{i=1}^{M}p(B_{i})\right)-1$ (22)

We exactly calculate simulated $\mathcal{F}_{\text{XEB}}(N_{q},d,s)$ using Eq.
(21) for each state $|\psi_{s}(N_{q},d)\rangle$, and then average over
$N_{c}50$ circuit instances with the same circuit depth to get
$\mathcal{F}_{\text{XEB}}(N_{q},d)$ in main text Fig. 3(c). The error bar in
the plot is the standard error over multiple circuit instances. For
experimental XEB (see main text Fig. 3(e)), we estimate
$\mathcal{F}_{\text{XEB}}(N_{q},d,s)$ using Eq. (22) for each state
$|\psi_{s}(N_{q},d)\rangle$, and then average over 50 circuit instances to get
estimated $\mathcal{F}_{\text{XEB}}(N_{q},d)$.

###  D.2 Data preprocessing

Out of 50 circuit instances for each depth, we randomly choose 35 circuit
instances as training circuits and the other 15 as testing circuits. For each
$|\psi_{s}(N_{q},d)\rangle$, we partition the $M$ measurement snapshots
sampled from the state into sets of set size $N$. For $N_{q}<36$ and the set
size of $N=10,000$, we have $35\times M_{s}/N=1750$ training sets per depth
$d$, $15\times M_{s}/N=750$ sets for testing each. For $N_{q}=36$, we have
$7000$ training sets and $3000$ testing sets. We combine all sets from the 35
(15) different circuit instances with the same circuit depth $d$ and label
them as the same class when generating the training (testing) dataset. We
train our models on data from two depths $d$ and $20$, labeled by $\hat{y}=0$
and $\hat{y}=1$ correspondingly.

As a baseline, we also train the model with shallow and deep depth both at
$d=20$. In this case, we expect the model to fail the classification task and
have accuracy at $50\%$. We randomly choose 35 circuit instances as a training
circuit from $d=20$, partition them into sets, and assign half of the sets to
shallow depth (label $\hat{y}=0$) and the other half to deep depth (label
$\hat{y}=1$). Testing datasets are constructed from the remaining 15 circuit
instances in the same way.

###  D.3 Training and testing procedure

We use PyTorch to train the models as binary classifiers to distinguish
shallow-depth circuit measurement outcomes from deep-depth circuit measurement
outcomes. We set the reference deep to be $d=20$ and vary the shallow depth.
We use Adam optimization with binary cross-entropy loss and also utilize the
dataset shuffling method, described in section C.3. Detailed model and
training hyperparameters are presented in Table 5 and 4. We employ a step
learning rate scheduler and Xavier normal initialization to enhance learning
further. We use the MSSAB module with $N_{s}=5$ to deal with large set size
$N=10,000$. We keep the model with the highest test accuracy during each
training run. The training curves are shown in SFig. 12.

Model hyperparameter
---
Model | QuAN50
Set size ($N$) | 10000
Number of mini-sets ($N_{\text{s}}$) | 5
Number of MSSAB layer ($L$) | 1
Number of $1\times 1$ Conv. channel ($n_{c}$) | $16$
Attention block: Hidden spatial dimension ($d_{h}$) | $16$
Attention block: Number of heads ($n_{h}$) | $4$
Attention block: activation function for residual connection | ReLU
Training hyperparameter
Optimizer |  `Adam`$(\beta_{1}=0.9,\;\beta_{2}=0.999,\;\epsilon=1\times 10^{-8})$
L2 coefficient | $5\times 10^{-5}$
Learning rate | $3.5\times 10^{-5}$
LR schedule |  `StepLR`$(\text{stepsize}=100\sim 200,\gamma=0.65)$
Epoch | 400
Dataset shuffle period | $10$
Batchsize | $20$
Initialization | `xavier_normal`
GPU | A100 (80GB)
Table 3: Model setting and hyperparameters used to train the model with random
quantum circuit data.

For the results presented in the main text, we perform $8$ independent
training runs on noiseless simulated data for each depth pair $(d,20)$, where
$d\in\\{4,6,\cdots,20\\}$. Therefore, for each $(d,20)$, we have 8 models. We
then test the trained models with simulated testing data from the same depth
pair $d,20$. Note that this is different from the hard-core Bose-Hubbard model
(section C.3) and toric code (section E.3) studies, where we validate models
trained on specific training parameter points across the entire phase space.
We then calculate the mean and the standard error of the mean of the accuracy
from the $8$ independent models for each $(d,20)$, as presented in the main
text Fig. 3(d). We then do testing on experimental data generated with
Google’s Sycamore processor. For each $(d,20)$, we apply the $8$ models
trained on simulated data to classify experimental data at the same depth pair
(see main text Fig. 3(g)).

| MLP | CNN | Transf.
---|---|---|---
| Conv(1, 16, 2, 1, BatchNorm) | Conv3d(1, 4, (500,2,2), (50,1,1), BatchNorm) | PE2d($d=16)$
Encoder | SLP(256, 16, Sigmoid) | Conv3d(4, 4, (50,2,2), (5,1,1), BatchNorm) | SAB(32, 16, 4)
| SLP(16, 16, Sigmoid) | Conv3d(4, 16, (5,2,2), (1,1,1), BatchNorm) |
Decoder | SLP(16, 16, Sigmoid) | SLP(1600, 1, Sigmoid) | SLP(400, 1, log softmax)
| SLP(16, 1, Sigmoid) |  |
| SLP*($N$, 1, Sigmoid) |  |
| SMLP | PAB | QuAN
| Conv(1, 16, 2, 1, BatchNorm) | Conv(1, 16, 2, 1, BatchNorm) | Conv(1, 16, 2, 1, BatchNorm)
Encoder | SLP(256, 48, Sigmoid) | SLP(256, 48, Sigmoid) | MSSAB*(256, 16, 4, $N_{s}$)
| SLP(48, 16, Sigmoid) | SLP(48, 16, Sigmoid) |
Decoder | SLP(16, 48, Sigmoid) | PAB*(16, 16, 4) | PAB*(16, 16, 4)
| SLP(48, 1) | SLP(16, 1, Sigmoid) | SLP(16, 1, Sigmoid)
| Sum*, Sigmoid |  |
Table 4:  Detailed model architectures used in main text Fig. 3(e). Data from
$N_{q}=25$ with $5\times 5$ geometry was used. Each layer’s arguments as
followed: Convolutional layer as Conv($n_{c,\text{in}}$, $n_{c,\text{out}}$,
kernel, stride, normalization). Single perception as SLP($d_{\text{in}}$,
$d_{\text{out}}$, activation). 2D Positional encoding as PE2d($d_{h}$).
SAB($d_{h,\text{in}}$, $d_{h,\text{out}}$, $n_{h}$) and
MSSAB($d_{h,\text{in}}$, $d_{h,\text{out}}$, $n_{h}$, $N_{s}$).
PAB($d_{h,\text{in}}$, $d_{h,\text{out}}$, $n_{h}$). * denotes the module
operates over a set dimension.

###  D.4 Machine learning details

Figure 12: (a) The learning curve of test classification accuracy compared at
different architectures using $N=10000$. (b) The learning curve of test loss
compared at different architectures. We use testing data from the system size
$N_{q}=25$ and shallow depth $d=8$. The solid line indicates the median of 5
independently trained models, and shaded regions show the minimum and maximum
test accuracy among 5 models at each epoch. Figure 13: Test classification
accuracy of QuAN50 with varying hyperparameters. (a) set size $N$, (b) number
of mini-sets $N_{s}$, (c) number of convolutional filters $n_{c}$, (d) size of
hidden dimension $d_{h}$, and (e) number of heads in multi-head attention
blocks $n_{h}$. Each gray dot represents each of the $5$ independently trained
models, and the purple star represents the maximal accuracy that we use in the
main text. The black solid lines with error bars represent the averaged
accuracy over $5$ models.

In SFig 12, we show the learning curves to check the model’s performance and
stability using test accuracy and loss curve as a function of epoch. Here, we
use data from depth pair $(8,20)$ for $N_{q}=25$. The test accuracy is
determined by the percentage of correctly classified test data that are from
either depth 8 or 20. Test loss is calculated through binary cross-entropy
between input true label $\hat{y}=0,1$ and model confidence $y(\mathbb{X})$
(see section A.1). The loss function quantifies how much machine confidence
deviates from the true label. The increasing test loss over epochs signals
overfitting to the training dataset and poor generalization to the testing
dataset. As previously discussed in the main text Fig. 3(e), only QuAN2 and
QuAN50 distinguish depth $8$ and $20$, while other architectures (MLP, CNN,
Transf., SMLP, PAB) have accuracies fluctuating near $50\%$. As the testing
loss curves for these models diverge or stay constant, this implies that the
model fails to distinguish between depth $8$ and $20$ data. Now, if we look at
the accuracy curves for QuAN2 and QuAN50, we observe both accuracy curves stay
higher than $50\%$ up to epoch$=400$ (see SFig. 12(a)). However, the testing
loss for QuAN2 is increasing, signaling overfitting of the training dataset,
while QuAN50 shows a decreasing trend in the loss function (see SFig. 12(b)).
Moreover, QuAN50 accuracy saturates to a higher value at an earlier epoch
compared to QuAN2 and also shows smaller variance across $5$ independent
training runs (shaded region). This demonstrates QuAN50 has good stability
with smaller dependence on randomness (model initialization, selection of
training and testing circuit index, random seeds). We conclude QuAN50 exhibits
stable and high performance. From now on, we exclusively employ QuAN50 for
random quantum circuit depth classification tasks.

In SFig. 13, we conduct a hyperparameter study for QuAN50 using the same data
(depth $8$ and $20$ with $N_{q}=25$ system, see Table 5). The hyperparameters
include set size $N$, number of mini-sets $N_{s}$, number of convolution
channel $n_{c}$, hidden dimension size $d_{h}$, and number of heads $n_{h}$
inside each attention block. As a baseline, we use
$(N,N_{s},n_{c},d_{h},n_{h})=(10000,5,16,16,4)$ and tune one type of
hyperparameter at a time. We record the highest testing accuracy for each
training. The plots show slight changes in average accuracy with varying
hyperparameters, except for the set size $N$. Unlike other hyperparameters,
set size $N=10000$ does not seem to have saturated the maximum accuracy yet.
However, further increasing the set size is impractical due to the
computational limit of the current GPU. Moreover, increasing the set size will
also result in a decrease in the total number of training sets (number of
training data points), which we encounter a risk of model overfitting to the
training data. Therefore, we conclude that our hyperparameter choice of
$(N,N_{s},n_{c},d_{h},n_{h})=(10000,5,16,16,4)$ yields optimal performance
given current computational resource.

Using the model QuAN50 with optimal hyperparameter setting, we train the model
using the data from different $(N_{q},d)$ to learn the evolution of state
complexity in random circuits as a function of depth $d$, shown in SFig. 14.
Our approach involves a multiple pairwise classification task, comparing data
from depth $d$ to that from depth $20$. The first classification task uses
simulated data without noise, whose result is shown in main text Fig. 3(d) and
SFig. 14(a). Next, we move on to classification using experimental data with
noise, shown in main text Fig. 3(g) and SFig. 14(b). To inspect experimental
data with QuAN, we utilize a model trained using depth $d$ and $20$ simulated
data, and then calculate test classification accuracy between depth $d$ and
$20$ experimental data. In the main text, we highlighted that the
classification accuracy for experimental data shows a similar trend to the
simulated data, exhibiting a sharp transition at depth $10$. However, an
exception occurs for system size $N_{q}=36$, where the accuracy trend differs
in distinguishing depth $4$ and $20$. To see whether a qualitative difference
exists for $(N_{q},d)=(36,4)$, we employ a new scheme of training and testing
both on the experimental data while keeping all other hyperparameter settings
the same (see SFig. 14(c)). We observe two key features from these curves.
Firstly, QuAN is readily distinguishing experimental data from depth $d$ and
depth $20$. Likely, this reflects the increasing degree of noise that comes
with circuit depth. Secondly, higher classification accuracy is seen on
$N_{q}=36$ data taken on a different day from the rest.

Figure 14: (a) Test classification accuracy by QuAN50 trained and tested with
both simulated data of $N_{q}$ qubits to distinguish depth $d$ and $20$. (b)
The accuracy of QuAN50 trained on simulated data and tested with experimental
data of depth $d$ and $20$. (c) The accuracy of QuAN50 trained and tested with
both experimental data.

##  Appendix E Toric code simulation

###  E.1 Data acquisition

We start with the $Z$-basis bit-string measurement from the deformed toric
code state with coherent noise that is now available in open source database
[10] as a part of the Ref. [11]. The database includes bit-strings obtained by
simulating measuring toric code deformed by coherent $X,Z$ noise with varying
strengths:

$\displaystyle|\psi(g_{X},g_{Z})\rangle=\frac{1}{\mathcal{N}}e^{-g_{X}\sum_{i}X_{i}-g_{Z}\sum_{i}Z_{i}}|\mathrm{TC}\rangle.$
(23)

The bit-strings available in the database are simulated and sampled using
projected entangled pair states (PEPS) on a $300\times 1000$ vertex square
lattice [11].

We then introduce incoherent noise through bit-flip with probability
$p_{\mathrm{flip}}$ (see main text Eq. (6)). Thus, the resulting bit-strings
are effectively sampled from mixed states $\rho(g_{X},g_{Z},p_{flip})$.

###  E.2 Data preprocessing

We construct training and testing datasets by transforming the $Z$-basis
measurements of the simulated mixed-state toric code into dual lattice sites
with dimensions of $300\times 1000$ $Z$-plaquette terms. For each state, we
extract $M_{s}=8134$ snapshots by slicing $300\times 1000$ dual lattice sites
into $6\times 6$ arrays, each containing 84 qubits. These snapshots represent
different quantum states, denoted as
$\\{\rho_{s}(g_{X},g_{Z},p_{\text{flip}})\\}$, where $s$ is index for
different state. The number of distinct states corresponding to a given
parameter set $(g_{X},g_{Z},p_{\text{flip}})$ ranges from 1 to 13. In our
analysis, we fix $g_{Z}=0.14$.

We create training and validation datasets at various bit-flip probability
$p_{\text{flip}}$, in both the topological phase
($p_{flip}\in\\{0,0.005,0.01,0.015,0.02\\}$) and trivial phase
($p_{\text{flip}}\in\\{0.3,0.305,0.31,0.315,0.32\\}$). To ensure an adequate
number of data points for training, we use different coherent noise points
($g_{Z}=0.14$ and $g_{X}\in\\{0,0.02,0.04,0.06,0.08\\}$) as topological or
trivial phase points (see hatched boxes in main text Fig. 4(c,d)).
Importantly, we strictly sample from the region where $g_{X}<g_{c}\approx
0.22$. In total, we have 200 distinct states, resulting in $M=N_{c}\times
M_{s}=200\times 8134$ snapshots per phase. Out of the 200 states, we randomly
allocate $150$ states (75%) for training and the remaining $50$ for
validation. For a given state $\rho_{s}(g_{X},p_{\text{flip}})$ with
$M_{s}=8134$ snapshots, we create a set $\mathbb{X}_{i}$ of set size $N$ by
partitioning snapshots into $\lfloor M_{s}/N\rfloor=127$ sets. Each snapshot
within this set is composed of the same state
$\\{\rho_{s}(g_{X},p_{\text{flip}})\\}$. For instance, we have
$150\times\lfloor M_{s}/N\rfloor=19050$ training sets, and $50\times\lfloor
M_{s}/N\rfloor=6350$ validation sets per phase for a set size of $N=64$.

After training the model, we test the model using testing dataset taken from
the entire phase space of coherent noise $0\leq g_{X}\leq 0.38$ and incoherent
noise $0\leq p_{\text{flip}}\leq 0.32$, which includes an intermediate region
(see SFig. 15). For the points inside the training region, we randomly choose
one state from the validation states. For the points outside of the training
region, we randomly choose one state from
$\\{\rho_{s}(g_{X},p_{\text{flip}})\\}$. Consequently, for the set size of
$N=64$, we use $\lfloor M_{s}/N\rfloor=127$ testing sets per each point
$(g_{X},p_{\text{flip}})$.

###  E.3 Training and testing procedure

Model hyperparameter
---
Model | SMLP, PAB, QuAN2
Set size ($N$) | 1$\sim$64
Number of mini-sets ($N_{s}$) | 1
Number of MSSAB layer ($L$) | 0,1
Conv. channel ($n_{c}$) | No convolution
Attention block: Hidden spatial dimension ($d_{h}$) | $16$
Attention block: Number of heads ($n_{h}$) | $4$
Attention block: activation function for residual connection | ReLU
Training hyperparameter
Optimizer |  `Adam`$(\beta_{1}=0.9,\;\beta_{2}=0.999,\;\epsilon=1\times 10^{-8})$
L2 coefficient | $5\times 10^{-5}$
Learning rate | $1\times 10^{-4}$
LR schedule |  `StepLR`$(\text{stepsize}=100,\gamma=0.65)$
Epoch | 200
Dataset shuffle period | $10$
Batchsize | $32768/N$
Initialization | `xavier_normal`
GPU | a100, 80GB
Table 5: Model setting and hyperparameters used to train the model for toric
code problem.

We use PyTorch to train and test the model to distinguish the topological and
trivial phases. We again minimize the binary cross entropy loss function using
Adam optimization, step learning rate scheduler, and dataset shuffling method
(see SM section C.3) to prevent the model from overfitting. Training
hyperparameters are listed in Table 3, and parameters of architectures are
listed in Table 6. We train three different architectures (ranging from SMLP
to QuAN2) and different set sizes (ranging from $N=1$ to $N=64$) independently
while keeping the remaining model parameters unchanged. For each architecture,
We also perform $10$ independent training to ensure stability. As listed in
Table 6, we no longer use the convolution layer that inspects spatial
correlation inside each snapshot. Instead, we make use of MLP as a function to
deal with nonlinear correlations between closed loops within the snapshot.

During training, we store the model with the highest validation accuracy and
later feed the testing dataset into the stored model. We obtain a phase
diagram from the machine confidence $y(\mathbb{X}_{i})$; we average the
machine output of confidence $y(\mathbb{X}_{i})$ over $127$ testing sets for
each point $(g_{X},p_{\text{flip}})$. We calculate the averaged confidence for
each of the 10 different models (which were trained independently) and obtain
the mean and the standard error for the averaged confidence.

| SMLP | PAB | QuAN2
---|---|---|---
Encoder | SLP(36, 48, Sigmoid) | SLP(36, 48, Sigmoid) | SLP(36, 48, Sigmoid)
| SLP(48, 16, Sigmoid) | SLP(48, 16, Sigmoid) | SLP(48, 16, Sigmoid)
|  |  | MSSAB*(36, 16, 4, 1)
Decoder | SLP(16, 48, Sigmoid) | PAB*(16, 16, 4) | PAB*(16, 16, 4)
| SLP(48, 1) | SLP(16, 1, Sigmoid) | SLP(16, 1, Sigmoid)
| Sum*, Sigmoid |  |
Table 6:  Detailed model architectures used in main text Fig. 2(e-g). Data
from $N_{q}=16$ with $4\times 4$ geometry was used. (Each layer’s arguments
are as follows: Convolution layer as Conv($n_{c,\text{in}}$,
$n_{c,\text{out}}$, kernel, stride, normalization). Single perception as
SLP($d_{\text{in}}$, $d_{\text{out}}$, activation). MSSAB($d_{h,\text{in}}$,
$d_{h,\text{out}}$, $n_{h}$, $N_{s}$). PAB($d_{h,\text{in}}$,
$d_{h,\text{out}}$, $n_{h}$)). * denotes the module operates on a set
dimension.

###  E.4 Benchmarking machine results

####  E.4.1 Benchmarking to locally error-corrected decoration (LED)

Figure 15:  (a,b) Phase diagrams of the toric code state with both coherent
noise (transverse field strength $g_{X}$ with fixed $g_{Z}=0.14$) and
incoherent noise $p_{\text{flip}}$ as a function of $Z$-loop tension $\alpha$
using (a) $1$ layer and (b) $3$ layers of locally error-corrected decoration
(LED) [11]. Yellow circles with vanishing dashed lines represent the known
thresholds at $p_{c}\approx 0.11$ and $g_{c}\approx 0.22$. (c,d) Phase diagram
constructed using averaged confidence $\bar{y}$ by QuAN2, presented in the
main text Fig. 4(c,d). The hatched regions mark the training data ranges.

Here, we present the extended results from QuAN trained with toric code data
through phase diagrams of the topologically non-trivial and trivial state,
which can be constructed from average confidence
$\bar{y}(g_{X},p_{\text{flip}})$ as a function of coherent noise $g_{X}$ and
incoherent bitflip noise $p_{\text{flip}}$ (see SFig. 15(c-d)). For
comparison, we present the phase diagram obtained using “locally error-
corrected decoration” (LED) (see SFig. 15(a-b), reproduced from Ref. [11]),
which classifies the topological phase based on the vanishing $Z$-loop tension
$\alpha$ after layers of operation. Here the loop tension $\alpha$ is defined
as $\langle Z_{\rm loop}(\gamma)\rangle=e^{-\alpha|\gamma|}$ for a loop
$\gamma$. As we discussed in the main text, we observe that QuAN2 sharpens the
transition and even saturates the known threshold $p_{c}\approx 0.11$ with
increasing set size $N$.

Figure 16: Sample complexity of (a) QuAN2 with varying set size $N$ and (b)
PAB and QuAN2 with set size $N=64$, using data from different
$p_{\text{flip}}$ along the $g_{X}=0$ axis. $N_{uc}=36$ is the number of unit
cells in snapshots. We use sample complexity of $\langle Z_{\rm
closed}\rangle$ as the baseline (marked in gray) [11]. Sample complexity
beyond the threshold (dashed line) is not defined.

We quantify sample complexity as another method to benchmark model performance
against LED, shown in SFig. 16. In our context, sample complexity is defined
as the number of samples (snapshots) required to confirm that the state is in
the topological phase with $95\%$ confidence. We aim to evaluate sample
complexity as a function of $p_{\text{flip}}$ along $g_{X}=0$ and see how the
required number of samples grows as we increase the incoherent noise level. We
utilize $13$ different states for each
$(g_{X},p_{\text{flip}})=(0,p_{\text{flip}})$ where $0.025\leq
p_{\text{flip}}\leq 0.3$, having a total $13\times 127=1651$ testing sets per
point. To ensure a fair comparison, we exclude points that overlap with the
training points.

To determine the sample complexity at a given point, we employ a t-test
between the model outputs from each set $y(\mathbb{X}_{i})$ and the
classification threshold at $y=0.5$. For each model, we first randomly sample
$D$ sets out of $1651$ sets from each $p_{\text{flip}}$, then feed in to
obtain $y(\mathbb{X}_{i})$, where $i$ runs from $1$ to $D$. (We randomly feed
each set into one stored model out of $10$ independently trained models with
equal probability.) We conduct the t-test with the null hypothesis that “those
$D$ sets are in trivial phase with an average confidence $\bar{y}\leq 0.5$”.
If the resulting p-value (probability of observing those outputs assuming
trivial phase) is less than $5\%$, we reject the null hypothesis, indicating
that the state is in the topological phase with over $95\%$ confidence. We
decrease the number of sets $D$, and $D^{*}$ is identified as the point where
the prediction fails to meet the $95\%$ confidence level. We define
$D^{*}\times N$ as the sample complexity, where $N$ is the size of each set.
This process is repeated $10$ times to ensure the stability of our
calculation, and the mean and standard error of $D^{*}\times N\times N_{uc}$
is plotted.

Upon obtaining sample complexities using different set sizes $N$ and
architectures (PAB, QuAN2), we compare them with the sample complexity using
bare Wilson loop $\langle Z_{\text{closed}}\rangle$ without LED [11] as a
baseline. Here, sample complexity refers to the number of samples required to
confirm $\langle Z_{\text{closed}}\rangle$ is non-zero with $95\%$ confidence,
calculated by $(2\sigma/\langle Z_{\text{closed}}\rangle)^{2}$ assuming
Gaussian distribution of loop expectation values. In SFig. 16, we present the
sample complexity at various $p_{\text{flip}}$ using models with different
hyperparameters. For QuAN2 with varying set sizes, we observe that the sample
complexity immediately increases exponentially for small set sizes (e.g.
$N=1$). However, with larger set sizes (e.g. $N=64$), sample complexity
remains relatively unchanged from its minimum value ($D^{*}\times N\times
N_{uc}=1\times 64\times 36=2304$) until $p_{\text{flip}}$ approaches phase
transition. Comparing this to the baseline (marked in grey in SFig. 16), the
advantage of using QuAN2 is clearly visible, especially with a large set size
($N=64$). Although the baseline sample complexity is lower for
$p_{\text{flip}}<0.035$, the sample complexity for QuAN2 with $N=64$ remains
constant even until $p_{\text{flip}}<0.09$. We hence conclude that QuAN2 is a
scalable method on a broader range of incoherent noise. The sample complexity
performance of PAB trends is similar to QuAN2 (see SFig. 16(b)), where both
QuAN2 and PAB show constant sample complexities on a broader range of
incoherent noise. This shows that PAB operation plays the main role in
maintaining a low sample complexity level.

####  E.4.2 Benchmarking to SMLP

Figure 17:  (a) Average confidence $\bar{y}$ by SMLP with varying set sizes
$N$, along the axis $g_{X}=0$ with different incoherent noise rates
$p_{\text{flip}}$. Green and pink stars represent the training points. (b)
Average confidence $\bar{y}$ by SMLP with varying set sizes $N$, along the
axis $p_{\text{flip}}=0$ with different coherent noise strength $g_{X}$. (c)
Average confidence $\bar{y}$ by PAB with varying set sizes $N$, along the axis
$g_{X}=0$ with different incoherent noise rates $p_{\text{flip}}$. (d) Average
confidence $\bar{y}$ by PAB with varying set sizes $N$, along the axis
$p_{\text{flip}}=0$ with different coherent noise strength $g_{X}$.

Now, we benchmark QuAN and PAB training to the simplest set-structured model,
SMLP (see SFig. 17). To see how the characterization of the topological phase
changes as we tune the architecture complexity and set size, we train three
different architectures (SMLP, PAB, QuAN2, see Table 6) with different set
sizes. In the main text Fig. 4(e-h), we made a comparison between QuAN2 of
various set sizes with SMLP($N=64$) and PAB($N=64$). Here, we would also like
to show set-size effects on SMLP and PAB, and make a comprehensive comparison
with QuAN. We draw the following two important conclusions from SFig. 17.
First, for both $g_{X}=0$ and $p_{\text{flip}}$ axis, we find that SMLP
prediction remains unchanged with increasing set size, and even introducing a
larger error bar due to a decrease in total training data points. In contrast,
PAB/QuAN sharpens the phase boundary with increasing set size (see also main
text Fig. 4(e-h) for QuAN results). This implies that even with a set
structure, treating every snapshot equally by averaging over snapshots does
not help to predict the topological phase. Moreover, comparing performances
between PAB and QuAN2, we conclude that the PAB module plays a critical role
in characterizing topological order. In the next section E.5, we will make a
comprehensive analysis of the PAB module.

###  E.5 Machine analysis: PAB as an importance-sampler

Figure 18: (a,c) Pooling attention score histogram from (a) topological state
with $(g_{X},p_{\text{flip}})=(0,0.05)$ (see main text Fig. 4(i)) and (c)
trivial state with $(g_{X},p_{\text{flip}})=(0,0.125)$. (b,d) Loop expectation
value $\langle Z_{\text{closed}}\rangle$ as a function of the loop perimeter
with high and low attention score in (b) topological
$(g_{X},p_{\text{flip}})=(0,0.05)$ (see main text Fig. 4(j)) and (d) trivial
state $(g_{X},p_{\text{flip}})=(0,0.125)$.

Here, we analyze the PAB module in more detail (see Fig. 4(i,j) in the main
text) and discuss the mechanism of the PAB in predicting the topological
phase. According to Eq. ((15)-(18)), the pooling attention output is given by

$\mathbb{p}_{\mu}=\sum_{\beta}s^{\prime\beta}\left(S_{\mu}+\sum_{\nu}V^{\prime}_{\mu\nu}\mathbb{z}_{\nu}^{\beta}\right)=\sum_{\beta}s^{\prime\beta}\mathbb{z}_{\mu}^{\prime\beta}\;\;\text{
with
}\;s^{\prime\beta}=\texttt{Softmax}\left[\sum_{\rho\lambda}\frac{1}{\sqrt{d_{h}}}\left({S}_{\rho}K^{\prime}_{\rho\lambda}\mathbb{z}_{\lambda}^{\beta}\right)\right],$
(24)

where $s^{\prime\beta}$ is the pooling attention score calculated from the
encoder output $\mathbb{z}^{\beta}$ with set index $\beta$.

Model training/testing hyperparameter
---
Model | PAB
Training set size ($N$) | $32$
Model index out of $10$ | $3$
Attention block: Head index | $3$
Testing data details
Testing set size ($N$) | $64$
Data points $(g_{X},p_{\text{flip}})$ | $(0,0.05)$ | $(0,0.125)$
Set number index ($i$) | $116$ | $43$
Table 7: Model hyperparameters and testing data details used in SFig. 18.

First, we demonstrate the attention score’s relation to $Z$-loop tension,
which serves as the order parameter of the topological phase. To this end, we
plot the histogram of pooling attention scores. Table 7 shows the detailed
hyperparameters used to obtain SFig. 18(a,c). After obtaining the pooling
attention score distribution for each set, we sample the highest $10$ ($\sim
15\%$) and lowest $10$ attention scores $s^{\prime\beta}$ and corresponding
snapshots $\\{B_{\beta}|s^{\prime\beta}\geq s_{\text{high}}\\}$ and
$\\{B_{\beta}|s^{\prime\beta}\leq s_{\text{low}}\\}$. For each snapshot
$B_{\beta}$ with $36$ $Z$-plaquette values, we calculate loop expectation
value $\langle Z_{\text{closed}}\rangle$ at different loop perimeter
$4L=4,8,12,16,20,24$ by multiplying $Z$-plaquettes inside the loop. SFig.
18(b,d) shows the mean and the standard error of the mean of loop expectation
value $\langle Z_{\text{closed}}\rangle$. The mean and standard error is taken
over $10$ high (low) attention score snapshots. For snapshots in the
topological phase, PAB assigns a high attention score on the snapshots with
vanishing loop tensions. The lower the loop tension is, the higher the
attention score is assigned to the snapshot. Since the module conducts
weighted sum over set index $\beta$ where weight is the attention score, PAB
acts as an automated importance sampler within a given set $\mathbb{X}_{i}$.

We then analyze how PAB makes a decision in testing. After Eq. (24), the
output is obtained through layer normalization, residual connection, and
single-layer perception:

$\displaystyle\mathbb{p}_{\mu}^{\prime}$
$\displaystyle=\texttt{LayerNorm}(\sum_{\beta}s^{\prime\beta}\mathbb{z}_{\mu}^{\prime\beta})$
(25) $\displaystyle y(\mathbb{X})$
$\displaystyle=\texttt{Sigmoid}\left(\sum_{\mu}W_{\mu}\;\texttt{LayerNorm}\left[\mathbb{p}_{\mu}^{\prime}+\texttt{rFF}_{\mu\nu}(\mathbb{p}^{\prime}_{\nu})\right]+b\right).$
(26)

where $\mathbb{z}=\texttt{Encoder}(\mathbb{x})$, and Encoder is a MLP for PAB
model. To simplify the process, we will compare between
$\mathbb{z}^{\prime\beta}_{\mu}$-vector and $W_{\mu}$-vector, ignoring latter
LayerNorm and residual layer which gives a minor shift in
$\mathbb{z}^{\prime\beta}_{\mu}$. SFig. 19 compares normalized encoder output
from high (low) attention score snapshots (from SFig. 18) with the final layer
weight matrix.

We notice that $\mathbb{z}^{\prime\beta}_{\mu}$-vector for high and attention
score snapshots are both parallel to the $W_{\mu}$-vector. Meanwhile,
$\mathbb{z}^{\prime\beta}_{\mu}$-vector from low attention score in trivial
phase is anti-parallel to the $W_{\mu}$-vector. This final Sigmoid activation
function then determines the output $y(\mathbb{X}_{i})$: if the inner product
between the importance-sampled encoder output $\mathbb{z}^{\prime\beta}_{\mu}$
and weight vector $W_{\mu}$ is positive (negative), it gives the final output
of topological phase $y>0.5$ (trivial phase $y<0.5$). When the snapshots with
high loop expectation values and high attention scores are no longer dominant,
the weighted average of encoder output becomes anti-parallel with the final
weight vector, and the machine prediction is no longer ‘topological’.

Figure 19: (a,b) Encoder output vector $\mathbb{z}$ for snapshots from high
(red) and low (blue) attention score in (a) topological
$(g_{X},p_{\text{flip}})=(0,0.05)$ and (b) trivial state
$(g_{X},p_{\text{flip}})=(0,0.125)$. (c) The final weight $W_{\mu}$ in final
perception before machine output (see Eq. (26)) shows a similar shape compared
to the encoder output vector with a high attention score. The dashed line
represents the value of bias $b$.

## References

  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, Attention is all you need, Advances in neural information processing systems 30 (2017).
  * Zaheer  _et al._ [2017] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Poczos, R. R. Salakhutdinov, and A. J. Smola, Deep sets, Advances in neural information processing systems 30 (2017).
  * Lee  _et al._ [2019] J. Lee, Y. Lee, J. Kim, A. Kosiorek, S. Choi, and Y. W. Teh, Set transformer: A framework for attention-based permutation-invariant neural networks, in  _International conference on machine learning_ (PMLR, 2019) pp. 3744–3753.
  * Karamlou  _et al._ [2024] A. H. Karamlou, I. T. Rosen, S. E. Muschinske, C. N. Barrett, A. Di Paolo, L. Ding, P. M. Harrington, M. Hays, R. Das, D. K. Kim, B. M. Niedzielski, M. Schuldt, K. Serniak, M. E. Schwartz, J. L. Yoder, S. Gustavsson, Y. Yanay, J. A. Grover, and W. D. Oliver, Probing entanglement in a 2D hard-core Bose–Hubbard lattice, Nature 10.1038/s41586-024-07325-z (2024).
  * Barrett  _et al._ [2023] C. N. Barrett, A. H. Karamlou, S. E. Muschinske, I. T. Rosen, J. Braumüller, R. Das, D. K. Kim, B. M. Niedzielski, M. Schuldt, K. Serniak, M. E. Schwartz, J. L. Yoder, T. P. Orlando, S. Gustavsson, J. A. Grover, and W. D. Oliver, Learning-based calibration of flux crosstalk in transmon qubit arrays, Phys. Rev. Appl. 20, 024070 (2023).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019), number: 7779 Publisher: Nature Publishing Group.
  * Hoke  _et al._ [2023] J. Hoke, M. Ippoliti, E. Rosenberg, D. Abanin, R. Acharya, T. Andersen, M. Ansmann, F. Arute, K. Arya, A. Asfaw, _et al._ , Measurement-induced entanglement and teleportation on a noisy quantum processor, Nature 622, 481 (2023).
  * Morvan  _et al._ [2023] A. Morvan, B. Villalonga, X. Mi, S. Mandra, A. Bengtsson, P. Klimov, Z. Chen, S. Hong, C. Erickson, I. Drozdov, _et al._ , Phase transition in random circuit sampling, arXiv preprint arXiv:2304.11119 (2023).
  * Developers [2023] C. Developers, Cirq (2023).
  * [10] N. Maskara, Enhancing detection of topological order by local error correction.
  * Cong  _et al._ [2024] I. Cong, N. Maskara, M. C. Tran, H. Pichler, G. Semeghini, S. F. Yelin, S. Choi, and M. D. Lukin, Enhancing detection of topological order by local error correction, Nature Communications 15, 1527 (2024).


# Condition for minimal Harmonic Oscillator Action

M. Moriconi<EMAIL_ADDRESS>Instituto de Física, Universidade Federal
Fluminense, Campus da Praia Vermelha, Niterói, 24210-340, RJ, Brazil.

###### Abstract

We provide an elementary proof that the action for the physical trajectory of
the one-dimensional harmonic oscillator is guaranteed to be a minimum if and
only if $\tau<\pi/\omega$, where $\tau$ is the elapsed time and $\omega$ is
the oscillator’s natural frequency.

##  I Introduction

Hamilton’s principle Fasano  establishes that the equations of motion for a
mechanical system can be obtained as a necessary condition for the action to
be locally stationary; that is, if the action is given by

$S=\int_{t_{1}}^{t_{2}}L(q,\dot{q},t)\,dt\,,$ (1)

where $L$ is the Lagrangian of the system ($q$ representing the generalized
coordinates), the condition $\delta S=0$ implies the equations of motion. We
are assuming that the system is holonomic and that the constraints have been
incorporated through a suitable choice of generalized coordinates.

The path that extremizes the action is often assumed to be the one that
minimizes it, but this is not guaranteed by the first-order variation. To
determine whether the action is a maximum or a minimum one must analyze the
second-order variation, similarly to what is done when studying critical
points of real functions of one variable. In fact, it can be quite tricky to
establish the minimality of the action for a given configuration. One of the
cases where this analysis can be performed rigorously is for the one-
dimensional harmonic oscillator (1D HO). A thorough analysis of this problem
is provided by Gray and Taylor.Gray

In order to understand the problem at hand, an analogy with a geometrical
problem is very useful. Consider the paths joining two points, $P$ and $Q$, on
the surface of a sphere of radius $R$. Without loss of generality, let us take
the point $P$ to be on the north pole. If $Q$ is close to $P$, we know that
the shortest path, as measured using the metric of the surface of the sphere,
is given by the arc of the great circle joining these two points, and that
this path is unique. As we move the point $Q$ farther away from $P$ this arc
will keep growing, until $Q$ is on the antipode to $P$, the south pole. Once
$Q$ is on the south pole, there is a surprising phenomenon—now there is an
infinite number of arcs joining $P$ and $Q$ with minimal length. For very good
reasons, therefore, the point $Q$ is called the kinetic focus. We should note
that the kinetic focus is more commonly known in the literature as the
conjugate point, and we refer the reader to the book by Gelfand and Fomin
Gelfand  for a readable treatment of this problem.

We see from this example that, given a point $P$ on the sphere, there exists a
point $Q$, its kinetic focus, such that the geodesics that are shorter than
the ones joining $P$ and $Q$ in fact minimize the path length. (The term
geodesic here refers to curves that satisfy the geodesic equation and are not,
necessarily, those with the shortest path length.) For geodesics that are
longer than those for the kinetic focus, it is no longer guaranteed that they
will minimize the path length. As a matter of fact, we know that an arc of a
great circle with length greater than $\pi R$ is not the shortest path, even
though it satisfies the geodesic equation.

In the action principle we are looking for paths that extremize the action,
and we can ask a similar question: given a solution of the equations of
motion—the true physical motion of a mechanical system—what are the conditions
that will guarantee that this solution minimizes the action? Notice that the
kind of boundary conditions that will be imposed here are crucial. We will
study the case where the space-time coordinates of the initial and final
points are given, which is necessary in writing the action principle.

This question was first analyzed by Legendre,Goldstine  who thought he had
solved it by looking at the second variation of the action $\delta^{2}S$. We
know that in order to guarantee that we have a minimum for the action, we must
have $\delta^{2}S>0$. Legendre’s conclusion for one-dimensional systems was
that, for a Lagrangian of the form $L(q,\dot{q},t)$, a necessary condition for
the second variation of the action to be non-negative should be
$\partial^{2}L/\partial\dot{q}^{2}\geq 0$.

Legendre’s condition is certainly necessary, for, intuitively, we can imagine
varying the solution of the equations of motion by adding a small perturbation
that changes the potential energy by a very small amount while changing the
kinetic energy by a large amount. If $\partial^{2}L/\partial\dot{q}^{2}$ can
be negative, this perturbation allows us to reduce the value of the action.
Only later, with the work of Jacobi,Goldstine  it was understood that one
needed to be careful about the existence of the kinetic focus.

For the harmonic oscillator there’s a similar phenomenon to the problem of
geodesics on a sphere. If we take the starting point to be $q_{0}=0$,
corresponding to the potential energy minimum, at $t_{0}=0$, and the end point
to be $q_{1}$, arbitrary, which is reached at time $\tau<\pi/\omega$, then
there’s only one solution for the equations of motion. On the other hand, if
we choose $\tau=\pi/\omega$, which corresponds to half of the period of
oscillation, we see that there is an infinite number of solutions connecting
$q_{0}=0$ at $t_{0}=0$ and $q_{1}=0$ at $t_{1}=\pi/\omega$, namely,
$q(t)=A\sin(\omega t)$ for any $A$. As shown by Gray and Taylor Gray , it is
possible to vary the path in such a way that the second variation of the
action vanishes. and and in the more general case where the starting and end
points are not at the minimum of the potential energy, as discussed by Gray
and Taylor Gray , to decide if we have a minimum or not.

We will show now that the solutions of the equation of motion of the 1D HO
connecting space-time points separated by a time interval less than
$\pi/\omega$, provides, in fact, a true minimum of the action. The Lagrangian
for the 1D HO is given by

$L=\frac{1}{2}m{\dot{q}}^{2}-\frac{1}{2}kq^{2}\,.$ (2)

Suppose we are given conditions $q(t_{i})$, $i=1,2$, and that we find a unique
solution for the equations of motion compatible with them, which we call
$q_{\star}(t)$. We should note that for such boundary conditions it is not
always the case that the solution is unique, as discussed by Gray.
ScholarpediaGray  Without loss of generality we take $t_{1}=0$ and
$t_{2}=\tau$. The question is, does this solution minimize the action for the
1D HO?

Normally, one approaches this question by writing an arbitrary function $q(t)$
as $q(t)=q_{\star}(t)+\eta(t)$, where $\eta(t_{i})=0$, $i=1,2$, and then
computing $S[q(t)]$ as

$\displaystyle S[q(t)]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t_{1}}^{t_{2}}\left[m(\dot{q}_{\star}+\dot{\eta})^{2}-k(q_{\star}+\eta)^{2}\right]dt$
(3) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t_{1}}^{t_{2}}\left(m{\dot{q}_{\star}}^{2}-kq_{\star}^{2}\right)dt+\frac{1}{2}\int_{t_{1}}^{t_{2}}\left(m{\dot{\eta}}^{2}-k\eta^{2}\right)dt$
$\displaystyle=$ $\displaystyle S[q_{\star}(t)]+S[\eta(t)]\,,$

where we have performed an integration by parts (the total derivative term
going to zero) and made use of the equation of motion. In order to show that
$S[q_{\star}(t)]$ is indeed the minimum action we must show that $S[\eta]\geq
0$ for non-zero $\eta(t)$, and equal to zero if and only if $\eta(t)=0$. The
traditional way to show thisGray  is to expand $\eta(t)$ in Fourier series and
then require $\delta^{2}S>0$. Due to the fact that
$\eta(t_{1})=\eta(t_{2})=0$, we can expand $\eta(t)$ in terms of $\sin(n\pi
t/\tau)$. The condition one obtains at the end is that $S[\eta(t)]\geq 0$ for
any $\eta(t)$, if and only if

$\tau<\frac{\pi}{\omega}\,,$ (4)

where $\omega^{2}=k/m$ is the oscillator’s natural frequency.

This approach offers a very powerful method, but it requires the use of heavy
mathematical machinery in order to establish the simple, but fairly important,
condition for $\tau$. Gray and Taylor Gray  present a derivation for the 1D
and 2D HO using the Fourier series method and a derivation using the general
theory of kinetic focus. In this note we establish Eq. (4) using only
elementary calculus.

Following a beautiful trick by Lax in a related variational problem,Lax  let
us write $\eta(t)=u(t)\sin(\pi t/\tau)$, with $u(t)$ arbitrary, since the sine
function guarantees that $\eta(t)$ vanishes for $t=0$ or $\tau$. Using this
function for $\eta(t)$ and defining $\phi=\pi t/\tau$, we find

$\displaystyle S[\eta(t)]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\tau}\left[m\left(\dot{u}\sin\phi+u\frac{\pi}{\tau}\cos\phi\right)^{2}-ku^{2}\sin^{2}\phi\right]dt$
(5) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\tau}\left[m\dot{u}^{2}\sin^{2}\phi+mu^{2}\left(\frac{\pi}{\tau}\right)^{2}\cos^{2}\phi\right.$
$\displaystyle\qquad\qquad\left.+2\frac{m\pi}{\tau}u\dot{u}\sin\phi\cos\phi-
ku^{2}\sin^{2}\phi\right]dt\,.$

We can then write

$2\frac{m\pi}{\tau}u\dot{u}\sin\phi\cos\phi=\frac{m\pi}{\tau}\left[\frac{d}{dt}\left(u^{2}\sin\phi\cos\phi\right)-\frac{\pi}{\tau}u^{2}\cos^{2}\phi+\frac{\pi}{\tau}u^{2}\sin^{2}\phi\right]\,,$
(6)

allowing us to integrate by parts, so that the action becomes

$S[\eta]=\frac{1}{2}\int_{0}^{\tau}\left[m\dot{u}^{2}\sin^{2}\left(\pi\frac{t}{\tau}\right)+m\left(\frac{\pi^{2}}{\tau^{2}}-\omega^{2}\right)u^{2}\sin^{2}\left(\pi\frac{t}{\tau}\right)\right]dt\,,$
(7)

where we have discarded the total derivative term, since it vanishes at $t=0$
and $t=\tau$.

The expression for the action given in Eq. (7) allows us to find the necessary
and sufficient conditions to guarantee that $S[\eta]$ is positive for all non-
zero $\eta$. If we choose $u(t)=\kappa$, a constant, we obtain

$S[\eta]=\frac{m}{2}\left(\frac{\pi^{2}}{\tau^{2}}-\omega^{2}\right)\kappa^{2}\int_{0}^{\tau}\sin^{2}\left(\pi\frac{t}{\tau}\right)dt=\frac{m}{2}\left(\frac{\pi^{2}}{\tau^{2}}-\omega^{2}\right)\kappa^{2}\frac{\tau}{2}\,,$
(8)

which, in order to be positive, requires

$\frac{\pi^{2}}{\tau^{2}}-\omega^{2}>0\quad\implies\quad\tau<\frac{\pi}{\omega}\,,$
(9)

establishing the fact that the Eq. (4) is a necessary condition. To establish
its sufficiency, we note that if the Eq. (4) is satisfied, then the integrand
in Eq. (7) is positive definite. Since the action is the sum of positive
definite terms, its minimum value is obtained when each term is $0$; that is,
for $u(t)=0$ and $S[q(t)]=S[q_{\star}(t)]$. The final result is that the
action is a minimum for time intervals less than half of the harmonic
oscillator period; for greater time intervals the action integral is a saddle
point, as defined by Gray and Taylor. Gray

###### Acknowledgements.

The author would like to thank the anonymous referees and the editor, whose
comments helped improve this paper considerably, and to Nivaldo Lemos, for a
critical reading.

## References

  * (1) A. Fasano and S. Marmi, Analytical Mechanics (Oxford University Press, 2006).
  * (2) C. G. Gray and E. F. Taylor, “When action is not least,” Am. J. Phys. 75, 434–458 (2007).
  * (3) I. M. Gelfand and S. V. Fomin, Calculus of Variations (Prentice-Hall Inc., 1963).
  * (4) See, for example, H. Goldstine, A History of the Calculus of Variations from the 17th through the 19th Century (Springer-Verlag, 1980), for the history of this problem and original references.
  * (5) Chris G. Gray, “Principle of least action,” Scholarpedia article; <http://www.scholarpedia.org/article/Principle_of_least_action>.
  * (6) P. Lax, “A short path to the shortest path,” American Mathematical Monthly 102, 158–159 (1995).