# Moving from continuous to discrete symmetry in the 2D XY model

Nouman Butt  Xiao-Yong Jin  James C. Osborn  Computational Science Division,
Argonne National Laboratory, 9700 S Cass Ave, Lemont IL 60439  Zain H. Saleem
Mathematics and Computer Science Division, Argonne National Laboratory, 9700 S
Cass Ave, Lemont IL 60439

###### Abstract

We study the effects of discretization on the U(1) symmetric XY model in two
dimensions using the Higher Order Tensor Renormalization Group (HOTRG)
approach. Regarding the $Z_{N}$ symmetric clock models as specific
discretizations of the XY model, we compare those discretizations to ones from
truncations of the tensor network formulation of the XY model based on a
character expansion, and focus on the differences in their phase structure at
low temperatures. We also divide the tensor network formulations into core and
interaction tensors and show that the core tensor has the dominant influence
on the phase structure. Lastly, we examine a perturbed form of the XY model
that continuously interpolates between the XY and clock models. We examine the
behavior of the additional phase transition caused by the perturbation as the
magnitude of perturbation is taken to zero. We find that this additional
transition has a non-zero critical temperature as the perturbation vanishes,
suggesting that even small perturbations can have a significant effect on the
phase structure of the theory.

##  I Introduction

While quantum computers have the potential to allow great advances in the
simulation of quantum systems, a major challenge will be making efficient use
of the limited resources available in near-term hardware. For systems with
continuous degrees of freedom, such as in common gauge theories, representing
the continuous field in a limited set of discrete states poses a challenge.
For quantum Hamiltonians, the mapping to an infinite set of discrete states
can be done through second quantization, but the set of states then needs to
be truncated to fit on a quantum computer. There have been several studies of
methods for discretizing and truncating field representations in lattice
models (both classical and quantum), including spin models [1, 2, 3, 4, 5, 6],
scalar fields [7, 8, 9, 10], and gauge theories [11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21] .

A similar truncation occurs when simulating these theories using tensor
networks [22]. Here we use this correspondence to study the effects of
truncation in two dimensional classical spin models as a proxy for mapping
models to a quantum computer. The example systems we will use are the $XY$
model and it’s discrete analog, the $Z_{N}$ symmetric clock model. The two-
dimensional XY model has a continuous $U(1)$ symmetry and exhibits a single
phase transition of the Berezinskii–Kosterlitz–Thouless (BKT) type [23, 24,
25]. The simplest discretization of the XY model, using $N$ discrete states,
produces the $Z_{N}$ clock-models, which break the $U(1)$ symmetry. This
discretization method has a major effect on the phase structure of the theory.
While the XY model is recovered in the $N\to\infty$ limit, one may need to go
to very large $N$ before the discretization effects become negligible at low
temperatures.

The $Z_{N}$ clock models have been studied extensively, and their phase
diagrams are known to be distinctly different from the XY model. For $N\leq
4$, they exhibit a single phase transition, though due to spontaneous symmetry
breaking, and not of the BKT type. For $N\geq 5$, there is clear evidence that
the models exhibit two phase transitions. The transition at higher temperature
becomes the BKT transition of the XY model as $N\to\infty$. For sufficiently
large, but finite $N$, it is also consistent with a BKT transition, even
though there is only a discrete symmetry. This is an example of an emergent
symmetry. An explanation for the emergent symmetry in the clock models was
provided based on a mapping of the classical model to a quantum Hamiltonian
[26, 27], which also suggested that the upper transition is BKT for all $N\geq
5$. However, numerical studies suggest that the $N=5$ transition may not be
the same type as for $N>5$, though may still be related [28, 29, 30]. The
nature of the lower temperature transition is also believed to be of BKT-type
based on central charge arguments [31, 32, 33]. While the upper critical
temperature remains fairly constant as $N\to\infty$, the lower critical
temperature moves towards zero in that limit [34]. It is the existence of the
lower critical temperature that largely distinguishes the $Z_{N}$ models for
$N\geq 5$ from the XY model.

Other discretization methods for the XY model are possible. In particular, the
formulation of the XY model as a tensor network, which can be done using a
character expansion, provides a convenient basis for truncation, and is
manifestly $U(1)$ invariant. Similar tensor network constructions using
character bases exist for gauge theories, and have been proposed to be useful
for mapping gauge theories to qubits while preserving gauge invariance [35].

The differences between truncation schemes for the XY model can be better
understood by examining the differences in their tensor network formulations.
It was shown that the tensor network formulation of the XY model preserves the
U(1) symmetry through the enforcement of selection rules [36, 37]. The
discrete version of the XY model, the $N$-state clock model, can be seen as a
version of the XY model where the infinite set of states are now folded into
the finite $N$ states in a periodic manner. The periodic nature of the $Z_{N}$
symmetry alters the selection rule and modifies the symmetry of the theory. To
further explore this we split the tensor construction into a core tensor,
centered on sites, and an interaction matrix, which connects neighboring
sites, and consider theories which mix the core and interaction from the
different models. From this we find that the core tensor, which enforces the
selection rule, is indeed the dominant factor on the phase structure of the
theory, as opposed to the interaction matrix.

We also examine the behavior of the phase structure when moving between the XY
and $Z_{N}$ models. As mentioned above, one way to transition from the clock
model to the XY model is to take $N$ large. An alternative path is to
introduce a symmetry breaking term into the XY model, that explicitly breaks
the symmetry down to $Z_{N}$. This is equivalent to the model studied
analytically by José, Kadanoff, Kirkpatrick, and Nelson (JKKN) [38, 39]. They
concluded that the $U(1)$ symmetry is unstable at low temperatures in response
to small values of the symmetry breaking field and the system develops a
second phase transition similar to the clock models. The perturbation has the
effect of smoothly transforming the XY core tensor into the $Z_{N}$ clock
model core tensor by the introduction of periodic terms. We are interested in
studying the emergence of the phase transition at low temperatures in the $XY$
model for small values of the symmetry breaking term. Our results suggest that
even small perturbations can have a significant effect on the low temperature
phase structure of the XY model. This may require caution when considering
simulations of truncated models to ensure that even small perturbations,
perhaps due to discretization effects or even simulation errors, don’t give
rise to new phases not present in the original theory.

In section II we present the tensor network representation of the $XY$ model,
its perturbed version and the $Z_{N}$ clock model, and elaborate on the
difference in the structure of the core tensor in each case. In section III we
compare numerical results for different truncations of $XY$ model based on the
character expansion and the $Z_{N}$ clock model. In section IV we show results
of mixing the core tensor and interaction matrix from different models, and
how it effects the phase structure. In section V we perform a detailed study
of the behavior of the lower temperature phase transition of the perturbed
model as the symmetry breaking field becomes smaller, to determine the effects
of small symmetry breaking on the phase diagram.

##  II Tensor network formulation of XY and $Z_{N}$ spin models

The Hamiltonians for the XY and $Z_{N}$ clock models can be written in the
common form

$H=-\sum_{\langle
ij\rangle}\cos(\theta_{i}-\theta_{j})-h\sum_{i}\cos(\theta_{i})$ (1)

where the first sum is over neighboring sites on a 2D periodic lattice and $h$
is the magnetic field. The angles, $\theta_{i}$, take continuous values in
$[0,2\pi)$ for the XY model and the discrete values $2\pi k/N$ with
$k\in{0,\ldots,N-1}$ for the $N$-state clock model.

The partition function for either model can be formulated as a tensor network
by separating the interaction term through an expansion in some basis, then
integrating (summing for $Z_{N}$) over the spin variables at each site [22,
40]. Details of this approach are given in Appendix A. The expansion produces
an interaction matrix on the links of the lattice while the integration
produces a core tensor on the sites of the lattice. The interaction matrix can
then be factored and absorbed into the core tensor to produce a single tensor
per site which describes the partition function of the system.

When defining the core tensor and interaction matrix, there is a choice of
where to put the local (single-site) terms in the Hamiltonian. One choice is
to keep the interaction matrix diagonal, which in the models considered here,
coincides with keeping only the interaction term
$H_{I}(\theta_{i},\theta_{j})=-\cos(\theta_{i}-\theta_{j})$ in the interaction
matrix. The choice of where to put the local terms has no effect on the
partition function or other physical results, and is only relevant when
studying the roles of the core tensor and interaction matrices on the phase
structure as done in section IV.

A common basis to use for the expansion of the interaction matrix is the
character basis, which, for the $U(1)$ symmetric XY model is just the phases
$\exp(in\theta)$. For the XY model, this gives the core tensor and interaction
matrix at inverse temperature $\beta$ of

$\displaystyle C^{XY}_{abcd}$ $\displaystyle=$ $\displaystyle
I_{a-b+c-d}(\beta h)$ (2) $\displaystyle M^{XY}_{ab}$ $\displaystyle=$
$\displaystyle I_{a}(\beta)\delta_{ab}$ (3)

where $I_{n}(x)$ is a modified Bessel function. The order of indices on the
core tensor, $C^{XY}_{abcd}$, is such that $a$ and $b$ (likewise $c$ and $d$)
correspond to opposite directions along the same dimension. Each of the
indices on the XY model core tensor extend from $-\infty$ to $\infty$. In
practice this needs to be truncated. Here we choose to truncate them
symmetrically in the range $-S,\ldots,S$ giving a total of $D=2S+1$ states for
each index. We will study the effect of this truncation for various values of
$D$ below.

For the $Z_{N}$ clock model, keeping only the interaction term
$H_{I}(\theta_{i},\theta_{j})$ in the interaction matrix, gives the core
tensor and interaction matrix

$\displaystyle C^{N}_{abcd}$ $\displaystyle=$
$\displaystyle\sum_{\ell}I_{a-b+c-d+\ell N}(\beta h)$ (4) $\displaystyle
M^{N}_{ab}$ $\displaystyle=$
$\displaystyle\sum_{k}I_{a+kN}(\beta)\delta_{ab}~{}.$ (5)

The indices on the clock model tensor are all finite, taking on $N$ distinct
values. Note that the difference in the XY and $Z_{N}$ tensors is that the
clock model has the same terms from the XY model, but they are folded around
the finite set of $N$ states. This is a consequence of the periodic nature of
the basis functions with discrete angles
$\exp(in\theta_{k})=\exp(in\theta_{k+N})$ for $\theta_{k}=2\pi k/N$ (see
Appendix A for details).

As detailed in Appendix A, the clock model tensors can also be constructed
using an infinite number of states which provides convenient approach to
directly compare with the XY model. In this case the tensors are

$\displaystyle C^{N\infty}_{abcd}$ $\displaystyle=$
$\displaystyle\sum_{\ell}I_{a-b+c-d+\ell N}(\beta h)$ (6) $\displaystyle
M^{N\infty}_{ab}$ $\displaystyle=$ $\displaystyle I_{a}(\beta)\delta_{ab}~{}.$
(7)

and the indices now extend from $-\infty$ to $\infty$. This is just the
interaction matrix from the XY model combined with the clock model core
tensor, extended to an infinite number of states on each index. This form of
the clock model provides a smooth interpolation between it and the XY model by
adjusting the terms in the core tensor, which is what the perturbed XY model
studied here will do.

For the perturbed $XY$ case (JKKN model) the following term is added to the
Hamiltonian

$\displaystyle\delta H=-h_{N}\sum_{i}\cos(N\theta_{i})$ (8)

which breaks the U(1) symmetry down to $Z_{N}$. This additional term modifies
the core tensor, giving

$\displaystyle C^{XYN}_{abcd}$ $\displaystyle=$
$\displaystyle\sum_{\ell}I_{a-b+c-d+N\ell}(\beta h)\frac{I_{\ell}(\beta
h_{N})}{I_{0}(\beta h_{N})}~{}.$ (9)

The normalization is chosen such that as $h_{N}\to\infty$ this reproduces the
$Z_{N}$ clock model core tensor, while at $h_{N}=0$ this is equivalent to the
XY model.

Given the interaction matrix and core tensor for a particular model, the site
tensor can be constructed as in Appendix A. By contracting a 2D lattice
network of site tensors, the partition function, $Z$, can be evaluated.
Observables can be obtained from derivatives of the partition function, either
numerically by finite differences of parameters, or by replacing a site tensor
with its derivative (an impurity tensor [41]). We use the impurity tensor
method to measure the magnetization. For the temperature derivatives (used in
the specific heat and cross-derivative) we use a the numerical derivative
obtained from a local polynomial fit.

The contraction of the tensors is performed numerically using the Higher Order
Tensor Renormalization Group (HOTRG) method [42]. This combines pairs of
neighboring site tensors to form a new blocked site tensor. As the size of a
tensor index grows above some limit, $D_{cut}$, we employ a Higher Order
Singular Values Decomposition (HOSVD) to truncate it back to $D_{cut}$ states.
This is distinct from the initial $D$ states in the truncation of the XY
model. In all cases we perform 20 blocking steps in each direction giving a
final volume of $V=2^{20}\times 2^{20}$ sites.

Figure 1: Specific heat, $C_{V}$, versus temperature, $T$, for the truncated
XY model with varying initial bond dimension $D$.

Figure 2: Specific heat, $C_{V}$, versus temperature, $T$, for the $Z_{N}$
clock model with varying $N$.

##  III Discrete approximations of continuous symmetry

Here we compare different ways to approximate the continuous XY model with a
fixed number of states in the tensor network representation. In particular we
compare the truncation of the XY model expanded in the character basis [22] to
the $Z_{N}$ clock model for the same number of states.

We first look at the effect of truncation on the specific heat,
$C_{V}=-T\partial_{T}^{2}F$, versus temperature, $T$, with varying number of
states. Here, $F=-(\ln Z)/(\beta V)$ is the free energy. The specific heat is
determined from the second derivative of a seven-point polynomial fit to the
free energy. The results in this section were obtained at $D_{cut}=40$.

In the case of the XY model, Fig. 2, for all values of the initial bond
dimension $D$, we see a peak in $C_{V}$ around $T\approx 1$ which is larger
than the BKT transition temperature of $T_{c}^{XY}\approx 0.89$ [43, 44, 45,
46]. While the peak of $C_{V}$ does not match the transition temperature (as
was previously known), it still happens to correlate well with the phase
structure of the theory (which we can see from other observables below). The
height of the $D=3$ peak is very different from the rest ($D\geq 5$). This is
likely related to the fact that the three state clock model is in a different
universality class than for $N\geq 5$. We also see $C_{V}$ rise again at low
temperature. This is more pronounced at small $D$ and the effect moves to
lower temperatures for larger $D$. At $D=15$, the specific heat matches that
from larger $D$ (up to $D=40$) very well down to $T=0.2$.

We can compare the results from the truncated XY model to those of the clock
models for the same number of states. In Fig. 2 we show the specific heat
versus temperature for the $N$-state clock models over the same range of $N$
as for $D$ before. Again we see that for three states, the behavior deviates
from the rest, but it is even more enhanced in this case. For $N\geq 5$ the
model develops a clear second peak in the specific heat at lower temperature.
The $N=5$ specific heat is closer to $N\geq 7$ than for $N=3$, but still
deviates fairly significantly. For $N\geq 7$ the peak around $T=1$ remains
consistent among the models with the main difference being the lower
temperature peak moving towards $T=0$ as $N$ increases. The deviations at
lower temperature are larger for the same number of states as for the
truncated XY case.

Figure 3: Magnetization per site, $\langle m\rangle$, versus temperature, $T$,
for the truncated XY model with magnetic field, $h=10^{-4}$, and varying
initial bond dimension $D$.

Figure 4: Magnetization per site, $\langle m\rangle$, versus temperature, $T$,
for the $Z_{N}$ clock model with magnetic field, $h=10^{-4}$, and varying $N$.

In Fig. 4 and Fig. 4 we compare the behavior of the average magnetization per
site, $\langle m\rangle=\partial_{h}F$, with a small magnetic field,
$h=10^{-4}$, between the two discretization schemes for the same number of
initial states. In the case of the truncated $XY$ model, the magnetization
never saturates to unity as we lower the temperature. However in the $Z_{N}$
clock models, we observe saturation at an intermediate temperature. As $N$
increases, the saturation to unity happens at lower temperatures.

Figure 5: Cross-derivative, $-\partial_{T}\langle m\rangle$, versus
temperature, $T$, for the truncated XY model with $h=10^{-4}$ and varying $D$.

.

Figure 6: Cross-derivative, $-\partial_{T}\langle m\rangle$, versus
temperature, $T$, for the $Z_{N}$ clock models with $h=10^{-4}$ and varying
$N$.

It is difficult to observe the lower temperature transition using magnetic
susceptibility. Instead we use a cross-derivative

$\displaystyle\frac{\partial^{2}F}{\partial h\partial
T}=-\frac{\partial\langle m\rangle}{\partial T}$ (10)

introduced in [33] to locate the lower temperature transition in the $Z_{N}$
models.

The cross-derivative of the magnetization for the XY model, Fig. 6, shows only
a single peak for all values of the initial number of states, as low as $D=3$.
Here the location of the $D=3$ peak differs from $D\geq 5$, in contrast to the
specific heat. This is comparable to results seen in the $S=1$ quantum XY
model [6]. For $D\geq 5$ the results are nearly identical for larger
temperatures and only show small changes in the value at lower temperatures
(shown in inset).

For the $Z_{N}$ models, the cross-derivative of the magnetization, Fig. 6,
clearly shows two peaks, indicating two phase transitions. The peak height for
the lower transition is smaller than the peak height for upper transition. As
$N$ is increased, we see that the lower peak moves towards lower temperatures.
In this case it is clear that the truncated $XY$ model using the character
basis approaches the XY model behavior quicker than the discrete $Z_{N}$
models for the same number of states.

Figure 7: Specific heat, $C_{V}$, versus temperature, $T$, for different
combinations of core tensor and interaction matrix from the truncated $XY$
model at $D=5$ and the $Z_{5}$ clock model. The combination with the XY
interaction matrix and $Z_{5}$ (periodic) core (squares) shows two peaks,
similar to the $Z_{5}$ model (triangles), while the $Z_{5}$ interaction and
$XY$ core (diamonds) shows no sign of a second phase transition.

Figure 8: Specific heat, $C_{V}$, versus temperature, $T$, for a mixed model
with the XY interaction matrix with a $Z_{5}$ core tensor generalized to
different numbers of initial states $D$. For $D=5$ this corresponds to the
same mixed model in Fig. 8 (squares). For large $D$ this approaches the
$Z_{5}$ clock model result (diamonds). The $XY$ model (circles) is also shown
for comparison.

##  IV Role of core tensor and interaction matrix

While both the core tensor and the interaction matrix differ between the XY
and clock models, it is expected that the core tensor, which implements the
selection rules responsible for charge conservation at $h=0$, is mainly
responsible for determining the symmetries of the theory [36]. We can check
this by combining the core tensor and interaction matrix from the different
theories and compare them to the original models.

In Fig. 8 we compare the different combinations of the core tensor and
interaction matrix from the $Z_{5}$ clock model and the XY model at $D=5$. The
combination with the XY interaction matrix and $Z_{5}$ (periodic) core tensor
clearly shows two peaks. This is similar to the $Z_{5}$ model, though the
lower temperature peak in the mixed case is much smaller. The mixed model with
a XY core and $Z_{5}$ interaction shows only a single peak within the
calculated temperature range, similar to the XY model. For these models it is
clear that the structure of the core tensor is the determining factor in the
overall phase structure of the theory.

In Fig. 8 we plot the specific heat of the $XY$ model with different numbers
of initial states and a suitably generalized $Z_{5}$ periodic core tensor. For
a large number of initial states the model reproduces the specific heat of
$Z_{5}$ clock model, and approaches the infinite dimensional representation of
the $Z_{N}$ models shown in eqs. (6) and (7). For $D=3$ there is only one
peak, while for $D\geq 5$ there are two peaks. The lack of a second peak in
the $D=3$ case can be understood due to there not being enough states to fully
feel the effects of the $Z_{5}$ periodicity in the core tensor. This
representation of the clock models provides a basis for studying the
transition between them and the XY model.

##  V Perturbed XY Model

In order to study the transition from the XY to clock models, we consider a
perturbed model by adding an additional term, eq. (8), which has an exact
$Z_{N}$ symmetry. This is equivalent to the model studied by José, Kadanoff,
Kirkpatrick, and Nelson (JKKN) [38, 39]. We choose a normalization, eq. (29),
for this model such that the $h_{N}\to\infty$ limit reproduces the $Z_{N}$
clock model. For numerical purposes, we use an approximation to the core
tensor in eq. (9). By writing it in the alternate form

$\displaystyle C^{XYN}_{abcd}$ $\displaystyle=$
$\displaystyle\sum_{k\ell}I_{k}(\beta h)\frac{I_{\ell}(\beta
h_{N})}{I_{0}(\beta h_{N})}\delta_{a-b+c-d+k+N\ell}~{},$ (11)

we can evaluate this for small $h$ by limiting the sum to $-2\leq k\leq 2$.
This form is accurate up to $O(h^{2})$, which is sufficient for our purposes.

For this model, we will always set $D=D_{cut}$, so that we are no longer
comparing the effects of truncation on the phase diagram. Instead we are
interested in looking at the effects of a symmetry breaking perturbation on
the continuous symmetry, separate from the effects of the truncation. One
could also combine the two and consider smaller truncation of the perturbed XY
model, but for simplicity, we do not consider that here. We note however, that
the truncated XY model at $D=5$ is already very close to the large $D$ limit
of the XY model (see e.g. Fig. 6), so we expect the main results here to carry
over even down to $D=5$.

We simulated this model at $N=5$, which is the smallest value of $N$ for which
a second phase transition appears in the clock model. We are interested in
seeing how the extra phase transition emerges for small $h_{5}$. The HOTRG
method used here requires an increasingly large $D_{cut}$ to get stable
results as we approach the continuous XY model, which limits the lower value
of $h_{5}$ which we could confidently simulate with our upper limit of
$D_{cut}=91$. In lieu of direct simulations at very small $h_{5}$ we instead
must extrapolate our results to small $h_{5}$.

Figure 9: Magnetization per site, $\langle m\rangle$, versus temperature, $T$,
for the perturbed XY model with magnetic field, $h=10^{-5}$, and varying
$h_{5}$.

Figure 10: Cross-derivative $-\partial_{T}\langle m\rangle$ versus
temperature, $T$, for the perturbed XY model with magnetic field, $h=10^{-5}$,
and varying $h_{5}$.

In Fig. 10 we plot the magnetization versus temperature with a small magnetic
field of $h=10^{-5}$ for different values of $h_{5}$ produced from HOTRG
simulations with $D_{cut}=40$. The magnetization profile smoothly interpolates
between the $XY$ and $Z_{5}$ models as the symmetry breaking field, $h_{5}$,
is varied in a continuous manner.

Again we will use the temperature derivative of the magnetization [33] to more
clearly identify the phase transitions. The temperature derivative is obtained
numerically from local polynomial fits to the magnetization. In Fig. 10 we
show $-\partial_{T}\langle m\rangle$ versus temperature for a range of
$h_{5}$, again with $h=10^{-5}$. Here we clearly see two peaks for all values
of $h_{5}$ shown. The upper peak occurs at a temperature, $T_{p2}$, which
shifts a little between $h_{5}=\infty$ (the $Z_{5}$ clock model) and
$h_{5}=1$, but is fairly stable below that. In contrast, the temperature of
the lower peak, $T_{p1}$, continues to move towards lower temperatures, and
also the peak height decreases. We will see below, however, that $T_{p1}$ does
not seem to go to zero as $h_{5}\to 0$, and is consistent with having a phase
transition at a temperature that is significantly away from zero as $h\to 0$,
for all values of $h_{5}$.

We obtain values for $T_{p1}$ and the corresponding peak height from a series
of fourth order polynomial fits to the magnetization versus temperature using
an interval of 11 data points, which represent a range in temperatures of
$\delta T=0.04$. This range was chosen as a good compromise between smoothing
out any small fluctuations in the data while preserving the sharpness of the
peaks as much as possible. Values of $T_{p1}$ are obtained from the zeros of
the second derivative of the fitted polynomials. Multiple values of $T_{p1}$
and the height are obtained by shifting the fit interval, and the values and
estimated errors are taken from the average and variance among the intervals
that contain the peak.

Figure 11: Lower peak temperature, $T_{p1}$, versus magnetic field, $h$, for
the perturbed XY model with varying $h_{5}$.

The results for $T_{p1}$ for a range of $h$ and $h_{5}$ values are shown in
Fig. 11. All these results were obtained using a bond dimension of $D=91$. As
noted above, the HOTRG results become less consistent as we approach smaller
values of $h$ and $h_{5}$. While we could get stable results for much smaller
$h_{5}$ within the range $10^{-7}\leq h\leq 10^{-5}$, we could only
confidently extract accurate peaks for $h_{5}$ down to 0.2. For $h_{5}<0.2$
the data becomes less smooth so that extracting accurate peaks becomes
difficult.

$h_{5}$ | $T_{c1}$ | $\gamma$
---|---|---
0.2 | 0.750(9) | 0.1257(12)
0.3 | 0.801(39) | 0.1195(19)
0.4 | 0.787(24) | 0.1336(20)
0.5 | 0.791(9) | 0.1295(65)
0.7 | 0.812(6) | 0.1315(28)
1.0 | 0.823(1) | 0.1353(25)
$\infty$ | 0.907(3) | 0.1220(51)
Figure 12: Table of fit results for the lower critical temperatures, $T_{c1}$,
from eq. (12), and peak height scaling exponent, $\gamma$, from eq. (15),
versus $h_{5}$.

For each value of $h_{5}$, we fit the data over the range $10^{-7}\leq h\leq
10^{-5}$ to the power-law scaling formula

$\displaystyle T_{c1}-T_{p1}\propto h^{\alpha}~{}.$ (12)

The fit values of the critical temperature, $T_{c1}$, are summarized in Fig.
12. Unfortunately, the error bars are too large to reliably extrapolate this
critical temperature to $h_{5}=0$.

Figure 13: Lower peak temperature, $T_{p1}$, versus $h_{5}$ for the perturbed
XY model with varying $h$.

An alternate approach is to first extrapolate $h_{5}\to 0$, then take the
$h\to 0$ limit. In Fig. 13 we use the same data from Fig. 11 to plot $T_{p1}$
versus $h_{5}$ for different values of $h$. We fit the data at fixed $h$ over
the range $0.2\leq h_{5}\leq 1.0$ to another power-law form

$\displaystyle T_{p1}-T^{*}_{p1}\propto h_{5}^{a}~{},$ (13)

where $T^{*}_{p1}$ is the fit result for the lower peak temperature
extrapolated $h_{5}\to 0$. These fits show good agreement across the whole
range in $h_{5}$.

Figure 14: Lower peak temperature, $T_{p1}$, extrapolated to $h_{5}=0$,
$T_{p1}^{*}$,versus magnetic field, $h$, for the perturbed XY model. The
smooth line is a fit to the power-law form in eq. (14).

The values of the extrapolated peak temperatures, $T_{p1}^{*}$, in the limit
$h_{5}\to 0$ are plotted in Fig. 14. For larger $h$, the extrapolated peak
temperature is consistent with zero. However for smaller $h$, the value of
$T_{p1}^{*}$ moves away from zero. A power-law fit to the form

$\displaystyle T^{*}_{c1}-T^{*}_{p1}\propto h^{\beta}~{},$ (14)

where $T^{*}_{c1}$ is the fit result for the lower peak temperature
extrapolated $h_{5}\to 0$ and $h\to 0$, is shown along with the data. The
extracted value of the critical temperature at $h=0$ and $h_{5}=0$ is
$T^{*}_{c1}=0.44(3)$. A similar result was obtained even when leaving out the
highest $h=4\times 10^{-6}$ point, so that point is not a significant factor
in constraining the curvature.

Figure 15: Lower peak height, $S_{p1}$, versus magnetic field, $h$, for the
perturbed XY model for varying $h_{5}$.

In order to establish whether the lower peak in the cross-derivative,
$-\partial_{T}\langle m\rangle$, corresponds to a true phase transition, we
look to see if the scaling of the peak height exhibits a critical behavior. In
Fig. 15 we plot the lower peak height, $S_{p1}$, versus $h$ for different
$h_{5}$. We see that the heights do increase as $h\to 0$. We fit the peak
height in the range $h\in[1e-7,4e-6]$ to a power-law scaling form

$\displaystyle S_{p1}\propto h^{-\gamma}~{}.$ (15)

The extracted exponents $\gamma$ are listed in Fig. 12. There doesn’t seem to
be any consistent change in $\gamma$ with $h_{5}$, and it appears that
$\gamma$ remains non-zero even as $h_{5}\to 0$. Thus the peaks appear to
correspond to a critical phase transition as $h\to 0$. This result is in
agreement with the RG analysis based conclusion from [38, 39], in which the
$h_{N}\to 0$ limit corresponds to a phase transition with non-zero critical
temperature.

##  VI Summary and Future Directions

We have studied the changes in phase structure when moving between the
continuous XY model and the discrete $Z_{N}$ clock model in two dimensions
using the Higher Order Tensor Renormalization Group approach. We compared
different ways to approximate the XY model with a fixed number of states in
the tensor network representation. We showed that the truncated character
expansion for the XY model approaches the continuous XY model behavior quicker
than than discrete $Z_{N}$ model with same number of states, as the number of
states increases.

We also examined the role of the core tensor and interaction matrix in
determining the phase structure of the XY and $Z_{N}$ models, and showed that
the core tensor plays the dominant role in the phase structure of these
models.

Further we explored a perturbed XY model (the JKKN model) which allows us to
interpolate between the XY and $Z_{N}$ models via the introduction of a
symmetry breaking term. This term breaks the $U(1)$ symmetry to $Z_{N}$ and
introduces a new phase transition (for $N\geq 5$) in addition to the well-know
BKT transition. We demonstrate that the cross-derivative peak corresponding to
this phase transition in the lower temperature region scales with volume, and
that the critical temperature doesn’t go to zero even in the limit of a
vanishingly small symmetry breaking field. This suggests that even small
symmetry-breaking perturbations can have a large effect on the phase structure
of theories. When discretizing models, and potentially when simulating them in
noisy environments, such as present quantum computers, one mush be careful to
avoid these possible effects.

While the perturbation considered here had an explicit $Z_{N}$ symmetry, which
may have had a dominant effect on the resulting phase structure, it would also
be interesting to explore the effect of perturbations with less symmetry to
see how the results are modified. For example, $N$ could be extended away from
integer values as was considered in the limit of a large perturbation here
[47, 48]. Determining the effect of small perturbations with varying degrees
of symmetry could be important to understanding the errors inherent in
simulations on resource limited and noisy near-term quantum simulators.

In this paper we have studied the XY model due to its simplicity and the fact
that its group-space discretizations, the clock models, exhibit notably
different phase structure. For future work it will be interesting to study the
discrete truncations of other continuous models such as the Heisenberg model,
O(N) vector model, the compact hyperbolic spin model [49], or gauge theories.
It would also be interesting to explore other expansion bases for
discretization to see their effects on the phase structure and if there are
alternatives that could produce better representations for a limited number of
states.

Acknowledgements N.B., X.-Y.J. and J.C.O. were supported through a QuantISED
grant by the U.S. Department of Energy, Office of Science, Office of High
Energy Physics. Z.S. is supported by the National Science Foundation under
Award No. 2037984 and DOE Q-NEXT. This research used resources of the Argonne
Leadership Computing Facility, which is a DOE Office of Science User Facility
supported under Contract DE-AC02-06CH11357.

##  Appendix A Tensor Network construction

Here we give a general treatment on the construction of the tensor network
formulation of spin model partition functions. This is a generalization of the
method used for the XY model [22].

We start with a Hamiltonian that is split into a nearest neighbor interaction
term, $H_{I}$, and local term, $H_{L}$, as

$\displaystyle H=\sum_{\langle
ij\rangle}H_{I}(\theta_{i},\theta_{j})+\sum_{i}H_{L}(\theta_{i})$ (16)

where the $\theta_{i}$ are variables on a site of the lattice that will be
integrated over, and $\langle ij\rangle$ represents all pairs of neighboring
sites.

Note that we can move local terms to or from the interaction term by defining

$\displaystyle H_{I}^{\prime}(\theta_{i},\theta_{j})$ $\displaystyle=$
$\displaystyle
H_{I}(\theta_{i},\theta_{j})+\frac{1}{2d}\left[F(\theta_{i})+F(\theta_{j})\right]$
(17) $\displaystyle H_{L}^{\prime}(\theta)$ $\displaystyle=$ $\displaystyle
H_{L}(\theta)-F(\theta)$ (18)

for arbitrary $F$. Setting $F=H_{L}$ would then then eliminate the local term
and put the whole model into the interaction term. This will be a convenient
form when discussing the $Z_{N}$ model below. For the XY model we will
consider the original form with the simpler interaction term. When
constructing the tensor network using a complete basis set, then choice of
convention for the local and interaction term does not matter. However, when
truncating the basis set, the different forms can give different results.

In constructing the tensor network formulation, the main step is to separate
the interacting variables by expanding the interaction term in some set of
basis functions, $f_{a}(\theta)$ and $g_{a}(\theta)$, as

$\displaystyle\mathrm{e}^{-\beta
H_{I}(\theta_{i},\theta_{j})}\approx\sum_{ab}M_{ab}f_{a}(\theta_{i})g_{b}(\theta_{j})$
(19)

with $M_{ab}$ a matrix characterizing the interaction. After collecting terms
with common integration variables we are left with a tensor on each site. On a
two dimensional square lattice this is given by

$\displaystyle C_{abcd}=\int d\theta~{}\mathrm{e}^{-\beta
H_{L}(\theta)}f_{a}(\theta)g_{b}(\theta)f_{c}(\theta)g_{d}(\theta)$ (20)

where $a,b$ and $c,d$ are pairs of opposing directions on the lattice.

The partition function is then a contraction of neighboring core tensors, $C$,
with the interaction matrices, $M$, placed in between. For convenience, one
typically defines a new tensor which incorporates the interactions into the
core to make evaluations simpler. One possibility is to take the SVD of the
interaction matrix as $M=U\Lambda V$ and then combine the factors
$\tilde{U}=U\sqrt{\Lambda}$ and $\tilde{V}=\sqrt{\Lambda}V$ with $C$ to get

$\displaystyle
T_{abcd}=C_{ijkl}\tilde{U}_{ia}\tilde{V}_{bj}\tilde{U}_{kc}\tilde{V}_{dl}~{}.$
(21)

Note that this can be viewed a performing a change of basis functions to make
$M$ the identity. This form is not unique though due to the gauge freedom in
tensor networks. Performing a similarity transformation on opposing pairs of
indices of the tensor still produces the same partition function.

###  A.1 XY model tensors

A standard basis used for a $U(1)$ spin variables, related to the character
expansion, is

$\displaystyle f_{a}(\theta)$ $\displaystyle=$
$\displaystyle\mathrm{e}^{ia\theta}$ (22) $\displaystyle g_{a}(\theta)$
$\displaystyle=$ $\displaystyle\mathrm{e}^{-ia\theta}~{}.$ (23)

If the interaction Hamiltonian is a function of $\theta_{i}-\theta_{j}$, then
this basis makes $M$ diagonal. For the standard XY model with
$H_{I}(\theta_{i},\theta_{j})=-\cos(\theta_{i}-\theta_{j})$, the interaction
matrix in the Fourier basis above, eqs. (22) and (23), is

$\displaystyle M^{XY}_{ab}=\delta_{ab}I_{a}(\beta)~{}.$ (24)

The core tensor for the XY model with a magnetic field and perturbation term

$\displaystyle H_{L}(\theta)=-h\cos(\theta)-h_{N}\cos(N\theta)$ (25)

is given by

$\displaystyle C^{XYN}_{abcd}$ $\displaystyle=$ $\displaystyle\mathcal{N}\int
d\theta\,\mathrm{e}^{\beta h\cos(\theta)+\beta
h_{N}\cos(N\theta)+i(a-b+c-d)\theta}~{}~{}~{}~{}~{}~{}~{}~{}$ (26)
$\displaystyle=$ $\displaystyle 2\pi\mathcal{N}\sum_{k\ell}I_{k}(\beta
h)I_{\ell}(\beta h_{N})\delta_{a-b+c-d+k+N\ell}$ (27) $\displaystyle=$
$\displaystyle 2\pi\mathcal{N}\sum_{\ell}I_{a-b+c-d+N\ell}(\beta
h)I_{\ell}(\beta h_{N})~{}.$ (28)

In order to match the $Z_{N}$ model (given below) in the large $h_{N}$ limit,
we have included a normalization factor

$\displaystyle\mathcal{N}=\frac{1}{2\pi I_{0}(\beta h_{N})}~{}.$ (29)

This scales the full partition function by a factor of $\mathcal{N}^{V}$.
Using the asymptotic form of the Bessel functions [50]

$\displaystyle I_{\ell}(\beta h_{N})\approx\frac{1}{\sqrt{2\pi\beta
h_{N}}}\mathrm{e}^{\beta h_{N}}~{},$ (30)

which is independent of $\ell$, we have

$\displaystyle\lim_{h_{N}\to\infty}C^{XYN}_{abcd}$ $\displaystyle=$
$\displaystyle\sum_{\ell}I_{a-b+c-d+N\ell}(\beta h)~{},$ (31)

which, as we will show below, is a form of the $Z_{N}$ clock model written
with an infinite basis set (as opposed to the typical finite $N$-state basis).

###  A.2 $Z_{N}$ clock model

There are several ways to construct the tensor network for the $Z_{N}$ clock
model. Here we present a few variations starting with $N$-state
representations in the discrete angle basis, and then in a character basis. We
also present an infinite state representation that coincides with the
$h_{N}\to\infty$ limit of the perturbed XY model. Lastly we consider it as an
approximation of the XY model using a specific basis to expand the interaction
term.

####  A.2.1 $N$-state representations

The simplest way to construct the $Z_{N}$ clock model tensor network is to
evaluate it in the basis of the discrete angles $\omega^{N}_{a}=2\pi a/N$.
This gives an interaction matrix of

$\displaystyle M_{ab}=\mathrm{e}^{\beta\cos(\omega^{N}_{a}-\omega^{N}_{b})}$
(32)

and a core tensor of

$\displaystyle C_{abcd}=\mathrm{e}^{\beta
h\cos(\omega^{N}_{a})}\delta_{ab}\delta_{ac}\delta_{ad}~{}.$ (33)

The interaction matrix is not diagonal in this basis. For comparison with the
XY model, it is convenient to make it diagonal by rotating to the character
(Fourier) basis giving

$\displaystyle M^{N}_{ab}$ $\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\sum_{k,\ell}\mathrm{e}^{i(a\omega^{N}_{k}-b\omega^{N}_{\ell})+\beta\cos(\omega^{N}_{k}-\omega^{N}_{\ell})}$
(34) $\displaystyle=$ $\displaystyle\delta_{ab}\sum_{k}I_{a+kN}(\beta)~{}.$
(35)

The core tensor in this basis becomes

$\displaystyle C^{N}_{abcd}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n}\mathrm{e}^{i(a-b+c-d)\omega^{N}_{n}+\beta
h\cos(\omega^{N}_{n})}$ (36) $\displaystyle=$
$\displaystyle\sum_{\ell}I_{a-b+c-d+N\ell}(\beta h)~{}.$ (37)

The normalizations in the above expressions were chosen for convenience.

####  A.2.2 Infinite state representation

One can get an infinite representation of the $Z_{N}$ clock model by using the
same basis, eqs. (22) and (23), and expansion of the interaction term as that
used in the XY model. This gives the same interaction matrix, eq. (24). The
core tensor is different than the XY model due to the discrete sum over
angles. Instead it takes the same form as the $Z_{N}$ core tensor in eq. (37),
except that all the external indices, $a$, $b$, $c$, and $d$, extend over the
range $-\infty\ldots\infty$ instead of just over $N$ states.

The $N$-state $Z_{N}$ interaction matrix contains the same set of Bessel
functions as the infinite state one, just with them folded over on to $N$
states. This is a consequence of the periodicity of the core tensor
$C^{N}_{a+N,b,c,d}=C^{N}_{a,b,c,d}$, and similarly for other indices.

####  A.2.3 Representation as approximation of XY model

One can also view the $Z_{N}$ model as an approximation of the XY model using
the set of $N$ basis functions

$\displaystyle f_{a}(\theta)=g_{a}(\theta)=\chi^{N}_{a}(\theta)$ (38)

with the rectangular step functions, $\chi^{N}_{a}(\theta)$, being 1 for
$-\omega^{N}_{1/2}<\theta-\omega^{N}_{a}\leq\omega^{N}_{1/2}\pmod{2\pi}$ and 0
otherwise. This basis divides the angles into $N$ disjoint parts, with each
basis function being constant over one part and zero on the others.

In this case we need to consider the alternate form of the interaction term,
$H_{I}^{\prime}$ from eq. (17) with $F=H_{L}$, so that the local term is moved
into the interaction. The expansion of the interaction term, eq. (19), can be
approximated using

$\displaystyle
M_{ab}=\mathrm{e}^{\beta\cos(\omega^{N}_{a}-\omega^{N}_{b})+(\beta
h/4)\left[\cos(\omega^{N}_{a})+\cos(\omega^{N}_{b})\right]}~{}.~{}~{}~{}~{}$
(39)

The core matrix in this basis is simply a copy tensor
$C_{abcd}=\delta_{ab}\delta_{ac}\delta_{ad}$. This form is similar to the
original form given for the clock model in eqs. (32) and (33), and the
interaction matrix can be diagonalized using the same rotation to the
character basis.

This basis provides a simple way to get a measure of the amount of error in
the truncation by examining the error in the approximation of the interaction
term. One could compare the error in the expansion of the interaction term
among different finite basis choices. This could then be used to choose a
basis with the smallest error.

## References

  * Milsted [2016] A. Milsted, Matrix product states and the non-abelian rotor model, Phys. Rev. D 93, 085012 (2016).
  * Alexandru  _et al._ [2019a] A. Alexandru, P. F. Bedaque, H. Lamm, and S. Lawrence (NuQS Collaboration), $\sigma$ models on quantum computers, Phys. Rev. Lett. 123, 090501 (2019a).
  * Singh and Chandrasekharan [2019] H. Singh and S. Chandrasekharan, Qubit regularization of the $o(3)$ sigma model, Phys. Rev. D 100, 054505 (2019).
  * Bruckmann  _et al._ [2019] F. Bruckmann, K. Jansen, and S. Kühn, O(3) nonlinear sigma model in $1+1$ dimensions with matrix product states, Phys. Rev. D 99, 074501 (2019).
  * Bhattacharya  _et al._ [2021] T. Bhattacharya, A. J. Buser, S. Chandrasekharan, R. Gupta, and H. Singh, Qubit regularization of asymptotic freedom, Phys. Rev. Lett. 126, 172001 (2021).
  * Zhang  _et al._ [2021] J. Zhang, Y. Meurice, and S.-W. Tsai, Truncation effects in the charge representation of the o(2) model, Phys. Rev. B 103, 245137 (2021).
  * Jordan  _et al._ [2012] S. P. Jordan, K. S. M. Lee, and J. Preskill, Quantum algorithms for quantum field theories, Science 336, 1130 (2012), https://www.science.org/doi/pdf/10.1126/science.1217069 .
  * Jordan  _et al._ [2014] S. P. Jordan, K. S. M. Lee, and J. Preskill, Quantum computation of scattering in scalar quantum field theories, Quantum Info. Comput. 14, 1014–1080 (2014).
  * Macridin  _et al._ [2018] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Digital quantum computation of fermion-boson interacting systems, Phys. Rev. A 98, 042312 (2018).
  * Klco and Savage [2019] N. Klco and M. J. Savage, Digitization of scalar fields for quantum computing, Phys. Rev. A 99, 052335 (2019).
  * Brower  _et al._ [2004] R. Brower, S. Chandrasekharan, S. Riederer, and U.-J. Wiese, D-theory: field quantization by dimensional reduction of discrete variables, Nuclear Physics B 693, 149 (2004).
  * Byrnes and Yamamoto [2006] T. Byrnes and Y. Yamamoto, Simulating lattice gauge theories on a quantum computer, Phys. Rev. A 73, 022328 (2006).
  * Zohar and Burrello [2015] E. Zohar and M. Burrello, Formulation of lattice gauge theories for quantum simulations, Phys. Rev. D 91, 054506 (2015).
  * Buyens  _et al._ [2017] B. Buyens, S. Montangero, J. Haegeman, F. Verstraete, and K. Van Acoleyen, Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks, Phys. Rev. D 95, 094509 (2017).
  * Bañuls  _et al._ [2017] M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Efficient basis formulation for ($1+1$)-dimensional su(2) lattice gauge theory: Spectral calculations with matrix product states, Phys. Rev. X 7, 041046 (2017).
  * Hackett  _et al._ [2019] D. C. Hackett, K. Howe, C. Hughes, W. Jay, E. T. Neil, and J. N. Simone, Digitizing gauge fields: Lattice monte carlo results for future quantum computers, Phys. Rev. A 99, 062341 (2019).
  * Alexandru  _et al._ [2019b] A. Alexandru, P. F. Bedaque, S. Harmalkar, H. Lamm, S. Lawrence, and N. C. Warrington (NuQS Collaboration), Gluon field digitization for quantum computers, Phys. Rev. D 100, 114501 (2019b).
  * Ji  _et al._ [2020] Y. Ji, H. Lamm, and S. Zhu (NuQS Collaboration), Gluon field digitization via group space decimation for quantum computers, Phys. Rev. D 102, 114513 (2020).
  * Bender and Zohar [2020] J. Bender and E. Zohar, Gauge redundancy-free formulation of compact qed with dynamical matter for quantum and classical computations, Phys. Rev. D 102, 114517 (2020).
  * Haase  _et al._ [2021] J. F. Haase, L. Dellantonio, A. Celi, D. Paulson, A. Kan, K. Jansen, and C. A. Muschik, A resource efficient approach for quantum and classical simulations of gauge theories in particle physics, Quantum 5, 393 (2021).
  * Ciavarella  _et al._ [2021] A. Ciavarella, N. Klco, and M. J. Savage, Trailhead for quantum simulation of su(3) yang-mills lattice gauge theory in the local multiplet basis, Phys. Rev. D 103, 094501 (2021).
  * Liu  _et al._ [2013] Y. Liu, Y. Meurice, M. P. Qin, J. Unmuth-Yockey, T. Xiang, Z. Y. Xie, J. F. Yu, and H. Zou, Exact Blocking Formulas for Spin and Gauge Models, Phys. Rev. D 88, 056005 (2013), arXiv:1307.6543 [hep-lat] .
  * Berezinskiǐ [1971] V. L. Berezinskiǐ, Destruction of Long-range Order in One-dimensional and Two-dimensional Systems having a Continuous Symmetry Group I. Classical Systems, Soviet Journal of Experimental and Theoretical Physics 32, 493 (1971).
  * Berezinskiǐ [1972] V. L. Berezinskiǐ, Destruction of Long-range Order in One-dimensional and Two-dimensional Systems Possessing a Continuous Symmetry Group. II. Quantum Systems, Soviet Journal of Experimental and Theoretical Physics 34, 610 (1972).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Sun  _et al._ [2019] G. Sun, T. Vekua, E. Cobanera, and G. Ortiz, Phase transitions in the Zp and U(1) clock models, Phys. Rev. B 100, 094428 (2019), arXiv:1907.09082 [cond-mat.str-el] .
  * Ortiz  _et al._ [2012] G. Ortiz, E. Cobanera, and Z. Nussinov, Dualities and the phase diagram of the p-clock model, Nucl. Phys. B 854, 780 (2012), arXiv:1108.2276 [cond-mat.stat-mech] .
  * Baek and Minnhagen [2010] S. K. Baek and P. Minnhagen, Non-kosterlitz-thouless transitions for the $q$-state clock models, Phys. Rev. E 82, 031102 (2010).
  * Borisenko  _et al._ [2011] O. Borisenko, G. Cortese, R. Fiore, M. Gravina, and A. Papa, Numerical study of the phase transitions in the two-dimensional z(5) vector model, Phys. Rev. E 83, 041120 (2011).
  * Baek  _et al._ [2013] S. K. Baek, H. Mäkelä, P. Minnhagen, and B. J. Kim, Residual discrete symmetry of the five-state clock model, Phys. Rev. E 88, 012125 (2013).
  * Blöte  _et al._ [1986] H. W. J. Blöte, J. L. Cardy, and M. P. Nightingale, Conformal invariance, the central charge, and universal finite-size amplitudes at criticality, Phys. Rev. Lett. 56, 742 (1986).
  * Affleck [1986] I. Affleck, Universal term in the free energy at a critical point and the conformal anomaly, Phys. Rev. Lett. 56, 746 (1986).
  * Chen  _et al._ [2020] Y. Chen, K. Ji, Z. Y. Xie, and J. F. Yu, Cross derivative of the gibbs free energy: A universal and efficient method for phase transitions in classical spin models, Phys. Rev. B 101, 165123 (2020).
  * Chatterjee  _et al._ [2018] S. Chatterjee, S. Puri, and R. Paul, Ordering kinetics in the $q$-state clock model: Scaling properties and growth laws, Phys. Rev. E 98, 032109 (2018).
  * Meurice  _et al._ [2020] Y. Meurice, R. Sakai, and J. Unmuth-Yockey, Tensor lattice field theory with applications to the renormalization group and quantum computing, (2020), arXiv:2010.06539 [hep-lat] .
  * Meurice [2019] Y. Meurice, Examples of symmetry-preserving truncations in tensor field theory, Phys. Rev. D 100, 014506 (2019).
  * Meurice [2020] Y. Meurice, Discrete aspects of continuous symmetries in the tensorial formulation of abelian gauge theories, Phys. Rev. D 102, 014506 (2020).
  * José  _et al._ [1977] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model, Phys. Rev. B 16, 1217 (1977).
  * José  _et al._ [1978] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Erratum: Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model, Phys. Rev. B 17, 1477 (1978).
  * Chen  _et al._ [2017] J. Chen, H.-J. Liao, H.-D. Xie, X.-J. Han, R.-Z. Huang, S. Cheng, Z.-C. Wei, Z.-Y. Xie, and T. Xiang, Phase transition of the q-state clock model: Duality and tensor renormalization, Chinese Physics Letters 34, 050503 (2017).
  * Zhao  _et al._ [2010] H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai, and T. Xiang, Renormalization of tensor-network states, Phys. Rev. B 81, 174411 (2010).
  * Xie  _et al._ [2012] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Coarse-graining renormalization by higher-order singular value decomposition, Phys. Rev. B 86, 045139 (2012).
  * Hasenbusch [2005] M. Hasenbusch, The two-dimensional xy model at the transition temperature: a high-precision monte carlo study, Journal of Physics A: Mathematical and General 38, 5869 (2005).
  * Yu  _et al._ [2014] J. F. Yu, Z. Y. Xie, Y. Meurice, Y. Liu, A. Denbleyker, H. Zou, M. P. Qin, J. Chen, and T. Xiang, Tensor renormalization group study of classical $xy$ model on the square lattice, Phys. Rev. E 89, 013308 (2014).
  * Jha [2020] R. G. Jha, Critical analysis of two-dimensional classical xy model, Journal of Statistical Mechanics: Theory and Experiment 2020, 083203 (2020).
  * Ueda and Oshikawa [2021] A. Ueda and M. Oshikawa, Resolving the berezinskii-kosterlitz-thouless transition in the two-dimensional xy model with tensor-network-based level spectroscopy, Physical Review B 104, 165132 (2021).
  * Hostetler  _et al._ [2021a] L. Hostetler, J. Zhang, R. Sakai, J. Unmuth-Yockey, A. Bazavov, and Y. Meurice, Clock model interpolation and symmetry breaking in O(2) models, PoS LATTICE2021, 353 (2021a), arXiv:2110.05527 [hep-lat] .
  * Hostetler  _et al._ [2021b] L. Hostetler, J. Zhang, R. Sakai, J. Unmuth-Yockey, A. Bazavov, and Y. Meurice, Clock model interpolation and symmetry breaking in o(2) models, Phys. Rev. D 104, 054505 (2021b).
  * Gubser  _et al._ [2016] S. Gubser, Z. H. Saleem, S. S. Schoenholz, B. Stoica, and J. Stokes, Nonlinear sigma models with compact hyperbolic target spaces, Journal of High Energy Physics 2016, 1 (2016).
  * [50] DLMF, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/10.30.E4, Release 1.1.5 of 2022-03-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.


#  LiteEFG: An Efficient Python Library for Solving Extensive-form Games

Mingyang Liu 1, Gabriele Farina1, Asuman Ozdaglar 1
1 LIDS, EECS, Massachusetts Institute of Technology
1<EMAIL_ADDRESS>

###### Abstract

LiteEFG is an efficient library with easy-to-use Python bindings, which can
solve multiplayer extensive-form games (EFGs). LiteEFG enables the user to
express computation graphs in Python to define updates on the game tree
structure. The graph is then executed by the C++ backend, leading to
significant speedups compared to running the algorithm in Python. Moreover, in
LiteEFG, the user needs to only specify the computation graph of the update
rule in a decision node of the game, and LiteEFG will automatically distribute
the update rule to each decision node and handle the structure of the
imperfect-information game.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Tour: Implementation of Counterfactual Regret Minimization (CFR)
    1. 3.1 Basics of Computation Graph
    2. 3.2 Construction of Computation Graph
    3. 3.3 Visitation of Infosets
      1. 3.3.1 Static Graph
      2. 3.3.2 Dynamic Graph
    4. 3.4 Loading the Game
    5. 3.5 Training
    6. 3.6 Evaluation
    7. 3.7 Debug
    8. 3.8 Various Baselines
  4. 4 Benchmark
  5. 5 Conclusion

##  1 Introduction

The successes of reinforcement learning in solving various games, including Go
(Silver et al., 2016, 2017), Atari (Mnih et al., 2013), and Dota2 (Berner et
al., 2019), have increased interest in developing scalable approaches for
finding equilibrium strategies in _extensive-form games (EFGs)_. Compared to
Markov games where the game state is fully observable, a big challenge in EFG
is that decisions are made according to partial observation of the game state.
It further results in additional hardness in computing the expected utility of
taking an action, unlike Q-values in Markov games, since the expected utility
depends on the distribution of hidden information. For instance, in Texas
Hold’em, the utility of raising the bid depends not only on the private hands
of the decision maker but also on those of her opponents, which are unknown
when making the decision.

The recent interest in computation methods for solving EFGs also sparked
activity in developing libraries for executing algorithms over the game tree.
A popular library is OpenSpiel (Lanctot et al., 2019), which provides various
game environments with an easy-to-use Python API for researchers to test their
algorithms. However, the algorithms operating over OpenSpiel are usually very
slow since they are executed via Python. This motivates the need to devise a
library that incorporates both simple Python API and efficient backend
execution, similar to how TensorFlow (Abadi et al., 2016) and PyTorch (Steiner
et al., 2019) operate.

To address the challenge above, we propose LiteEFG, an open-source library
with a simple Python API and an efficient C++ backend for solving EFGs, which
is much faster than pure Python execution. With LiteEFG, researchers need to
define the update-rule at a decision node in a similar way as defining neural
networks with TensorFlow or Pytorch, then LiteEFG will automatically
distribute that update-rule to each individual decision node and handle the
relationship between different decision nodes automatically. Compared to
OpenSpiel (Lanctot et al., 2019), LiteEFG is simpler and faster when solving
tabular games (games that may fit into the computer memory). In experiments,
the classical baseline, Counterfactual Regret Minimization (Zinkevich et al.,
2007), implemented by LiteEFG is about $100\times$ faster than that of
OpenSpiel.

In LiteEFG, researchers only need to specify the computation graph of the
algorithm via Python. Then, the computation graph will be executed via C++,
which provides acceleration by several orders of magnitude. Moreover, due to
the imperfect information of EFG, the game states and decision nodes no longer
coincide, which complicates the implementation of the algorithm. To simplify
the issue, LiteEFG will automatically aggregate the information from different
game states belonging to the same decision node 111In EFGs, since the game is
partially observable, some game states may not be differentiable for the
decision maker. For instance, in Texas Hold’em, the game states that only
differ in the private hands of the opponents are not differentiable for the
decision maker. Then, we call those game states belonging to the same decision
node, because the decision made at those game states should be the same., so
users only need to specify the update-rule for the decision node, without
concerning with the aggregation process.

##  2 Preliminaries

In this section, we will introduce the preliminaries of EFGs. We use
$\Delta^{m}\coloneqq\left\\{\bm{x}\in[0,1]^{m}\colon\sum_{i=1}^{m}x_{i}=1\right\\}$
to denote the $m-1$ dimensional probability simplex. For a discrete set
$\mathcal{C}$, we use $|\mathcal{C}|$ to denote its cardinality. For any real
number $x\in\mathbb{R}$, $[x]^{+}\coloneqq
x\cdot\operatorname*{\mathds{1}}_{x\geq 0}$, which is $x$ when $x\geq 0$ and
$0$ otherwise.

##### Basics of extensive-form games.

In an $N$-player EFG, we use $[N]\coloneqq\left\\{1,2,...,N\right\\}$ to
denote the set of all players. Optionally, a fictitious player—called the
_chance player_ —is introduced to model stochastic events, such as a random
draw of cards or dice roll.

The game is a tree structure and we use $\mathcal{H}$ to denote the set of all
nodes in the game. For each $h\in\mathcal{H}$, one of the players among
$\left\\{c\right\\}\cup[N]$, where $c$ is the chance player, will take actions
at $h$. We use $p(h)$ to denote the player acting at $h$, and say that “$h$
belongs to $p(h)$”. Then, player $p(h)$ can choose an action in the action set
$\mathcal{A}_{h}$ and the process will be repeated until the game reaches a
terminal node $h^{\prime}\in\mathcal{Z}$, where $\mathcal{Z}$ is the set of
all terminal nodes. The utility for each player $i\in[N]$ is denoted by
$\mathcal{U}_{i}:\mathcal{Z}\to[0,1]$.

##### Information Set.

In an imperfect-information game, a player $i\in[N]$ may not be able to
distinguish all the nodes of the tree, that is,
$\left\\{h^{1},h^{2},...,h^{k}\colon\forall j=1,2,...,k,p(h^{j})=i\right\\}$.
For example, in the two-player Texas Hold’em, each player cannot distinguish
nodes that only differ because of her opponent’s private hand, since the hand
of the opponent is not revealed. To model imperfect information, a partition
of each player’s nodes—called the player’s _information partition_ —is
introduced. The elements of the partition, called _information set_ (or
infoset for short), denote nodes that are indistinguishable to the player when
acting at any of them. ${\mathcal{S}}_{i}$ is used to denote the set of all
infosets of player $i\in[N]$. For simplicity, for each node $h$, we use $s(h)$
to denote the infoset that $h$ belongs to and extend the player indicator $p$
to infosets, _i.e._ $p(s(h))\coloneqq p(h)$.

##### Strategy of Players.

The strategy of player $i\in[N]$ can be written as
$\pi_{i}(\cdot{\,|\,}s)\in\Delta^{|\mathcal{A}_{s}|}$ for each infoset
$s\in{\mathcal{S}}_{i}$, where $\mathcal{A}_{s}$ is the action set of infoset
$s$.222Since nodes in $s$ cannot be differentiated by $i$, the action set of
each node in $s$ must be the same. Therefore, infosets are also called
_decision nodes_ , since players make decisions conditioned on the infoset.

##### Reach Probability and Sequence-form Strategy.

For any strategy $\pi_{i}$ of player $i\in[N]$, we can define the reach
probability of player $i$ as $\mu_{i}^{\pi_{i}}(h^{1}\to h^{2})$ as the reach
probability from node $h^{1}$ to node $h^{2}$, to which we only count the
probability contributed by player $i$ while applying $\pi_{i}$. Moreover, we
can define $\mu_{i}^{\pi_{i}}(s\to h)$ as the reach probability from infoset
$s\in{\mathcal{S}}_{i}$ to node $h$, when only counting the probability
contributed by $\pi_{i}$. We use $\mu_{i}^{\pi_{i}}(\emptyset\to h)$
($\mu_{i}^{\pi_{i}}(\emptyset\to s)$) as the reach probability to $h$ ($s$)
from the root of the game. With reach probability, we slightly abuse the
notion $\mu$ to denote the sequence-form strategy
$\mu_{i}^{\pi_{i}}(s,a)\coloneqq\mu_{i}^{\pi_{i}}(\emptyset\to
s)\pi_{i}(a{\,|\,}s)$.

##  3 Tour: Implementation of Counterfactual Regret Minimization (CFR)

In this section, we will introduce how to use LiteEFG.

###  3.1 Basics of Computation Graph

LiteEFG is based on the computation graph, in which a vector is stored at each
node and users need to define the relationships between graph nodes. For
instance, node $A$ equals node $B$ plus node $C$. Then, every time the user
updates the graph, the variables at each graph node will be updated according
to the predefined relationship.

In LiteEFG, the user need to define the computation graph for an infoset
first. Then, the graph will be copied to each infoset in
$\bigcup_{i\in[N]}{\mathcal{S}}_{i}$. Therefore, all infosets share the same
relationship between graph nodes, while the variables stored in the graph node
of each infoset are independent.

In LiteEFG, the user can define a node by LiteEFG.function(...) with some
function of LiteEFG 333The full API list can be found in
https://github.com/liumy2010/LiteEFG., such as LiteEFG.sum and LiteEFG.exp. In
this case, LiteEFG will create a new node to store the outcome of the function
and return that node. Alternatively, the user can update the variable at a
node by x.inplace(LiteEFG.function(...)). In this case, LiteEFG will not
create a new node. Instead, the outcome of the function will be stored at node
x, and replace the original variable at x.

###  3.2 Construction of Computation Graph

In this section, we will introduce the construction of the computation graph
for Counterfactual Regret Minimization (CFR) (Zinkevich et al., 2007), one of
the most prominent algorithms for solving EFGs, with LiteEFG.

###  3.3 Visitation of Infosets

Figure 1: At the very beginning, when the computation graph for the
environment is determined by LiteEFG.Environment.set_graph, the static
backward nodes and static forward nodes will be executed sequentially. Then,
each time LiteEFG.Environment.update is called, the dynamic backward nodes and
dynamic forward nodes will be executed.

LiteEFG provides four types of graph nodes.

  * •

LiteEFG.backward(is_static=True): Static backward nodes. These nodes will be
executed at initialization (ahead of the execution of any other nodes). To
execute the static backward nodes, infosets will be visited in the reversed
breadth-first order and the corresponding static backward nodes will be
executed.

  * •

LiteEFG.forward(is_static=True): Static forward nodes. These nodes will be
executed at initialization (ahead of the execution of any dynamic nodes, but
after the static backward nodes). To execute the static backward nodes,
infosets will be visited in the breadth-first order and the corresponding
static forward nodes will be executed.

  * •

LiteEFG.backward(is_static=False): Dynamic backward nodes. These nodes will be
executed every time the function LiteEFG.Environment.update is called. To
execute the dynamic backward nodes, infosets will be visited in the reversed
breadth-first order and the corresponding dynamic backward nodes will be
executed.

  * •

LiteEFG.forward(is_static=False): Dynamic forward nodes. These nodes will be
executed every time the function LiteEFG.Environment.update is called. To
execute the dynamic backward nodes, infosets will be visited in the breadth-
first order and the corresponding dynamic forward nodes will be executed.

The order is also illustrated in Figure 1.

####  3.3.1 Static Graph

Figure 2: Implementation of Counterfactual Regret Minimization (CFR)
(Zinkevich et al., 2007) by LiteEFG. The green box displays the definition of
static variables that will only be updated once at initialization. The red box
displays the variables that will be updated every time updating the graph.

The update-rule of CFR is as follows. For any player $i\in[N]$, an infoset
$s\in{\mathcal{S}}_{i}$, and action $a\in\mathcal{A}_{s}$, the update-rule at
timestep $t\in\left\\{1,2,...,T\right\\}$ is

$\displaystyle\text{CF}^{(t+1)}_{i}(s,a)=\sum_{h\in\mathcal{Z}}\mathcal{U}_{i}(h)\frac{\mu_{i}^{\pi_{i}^{(t+1)}}(s\to
h)}{\pi_{i}^{(t+1)}(a{\,|\,}s)}\prod_{j\in\left\\{c\right\\}\cup[N]\colon
j\not=i}\mu_{j}^{\pi_{j}^{(t+1)}}(\emptyset\to h)$ $\displaystyle
R_{i}^{(t+1)}(s,a)=R_{i}^{(t)}(s,a)+\text{CF}_{i}^{(t+1)}(s,a)-\mathbb{E}_{a^{\prime}\sim\pi^{(t+1)}_{i}(\cdot{\,|\,}s)}\left[\text{CF}_{i}^{(t+1)}(s,a^{\prime})\right]$
(3.1)
$\displaystyle\pi^{(t+2)}_{i}(a{\,|\,}s)=\begin{cases}\frac{[R_{i}^{(t+1)}(s,a)]^{+}}{\sum_{a^{\prime}\in\mathcal{A}_{s}}[R_{i}^{(t+1)}(s,a^{\prime})]^{+}}&\sum_{a^{\prime}\in\mathcal{A}_{s}}[R_{i}^{(t+1)}(s,a^{\prime})]^{+}>0\\\
\frac{1}{|\mathcal{A}_{s}|}&\sum_{a^{\prime}\in\mathcal{A}_{s}}[R_{i}^{(t+1)}(s,a^{\prime})]^{+}=0\end{cases}$

From Equation 3.1, to update CFR, we need to maintain two variables in each
infoset $s$, the strategy $\bm{\pi}_{i}(\cdot{\,|\,}s)\in\Delta^{|A_{s}|}$,
and the regret buffer $R(s,\cdot)\in\mathbb{R}^{|A_{s}|}$. In Figure 2, the
green code block defines the static backward nodes of CFR algorithm. In line
7, we define expectation, which is the placeholder for
$\mathbb{E}_{a^{\prime}\sim\pi^{(t+1)}_{i}(\cdot{\,|\,}s)}\left[\text{CF}_{i}^{(t+1)}(s,a^{\prime})\right]$.
We will discuss the necessity of such placeholder in Section 3.3.2. Line 8 and
9 initialize $\bm{\pi}_{i}(\cdot{\,|\,}s)$ as uniform distribution over
$\Delta^{|\mathcal{A}_{s}|}$ and $R(s,\cdot)$ as zero vector individually. The
vector at self.action_set_size is a scalar equivalent to $|\mathcal{A}_{s}|$
for each infoset $s$.

####  3.3.2 Dynamic Graph

In this section, we will focus on the  red code block in Figure 2.

##### Line 15.

Line 15 displays the usage of LiteEFG.aggregate.

$\displaystyle{\tt LiteEFG.aggregate(}$ $\displaystyle{\tt
x,aggregator\\_name:["sum","mean","max","min"],}$ $\displaystyle{\tt
object\\_name:["children","parent"]="children",}$ $\displaystyle{\tt
player:["self","opponents"]="self",padding=0.0)}$

  * •

x: Indicates the graph node from which the function aggregates information.

  * •

aggregator_name: Method to aggregate the information.

  * •

object_name: Aggregate information from parent / children of the current
infoset. By default, it is “children".

  * •

player: Aggregate information from the infosets belong to self / opponents.
For instance, when object_name="children", for an infoset
$s\in{\mathcal{S}}_{i}$, we will aggregate the information from the children
$s^{\prime}$ of $s$, with $p(s^{\prime})=i$ if `player="self"` and
$s^{\prime}$ with $p(s^{\prime})\not=i$ otherwise.

  * •

padding: If the parent / children of the current infoset does not exist,
return padding.

When the object_name is “children”, for each action $a\in\mathcal{A}_{s}$ in
infoset $s$, there may be several subsequent children infoset. aggregate will
first concatenate the variable at x at all those infosets together. Then,
LiteEFG will call the aggregator function specified in aggregator_name to
aggregate that vector to a single scalar. Finally, LiteEFG will concatenate
the scalar at each action $a\in\mathcal{A}_{s}$ to a vector in
$\mathbb{R}^{|\mathcal{A}_{s}|}$.

When the object_name is “parent”, suppose $(s^{\prime},a^{\prime})$ is the
parent sequence of $s$. If the variable $\bm{v}$ stored at x at $s^{\prime}$
is in $\mathbb{R}^{|\mathcal{A}_{s^{\prime}}|}$, then $v_{a^{\prime}}$ will be
returned. Or if $\bm{v}$ will be returned if it is a scalar. Otherwise, an
error occurs.

##### Line 16.

Line 16 computes
$\mathbb{E}_{a^{\prime}\sim\pi^{(t+1)}_{i}(\cdot{\,|\,}s)}\left[\text{CF}_{i}^{(t+1)}(s,a^{\prime})\right]$
to replace the old value
$\mathbb{E}_{a^{\prime}\sim\pi^{(t)}_{i}(\cdot{\,|\,}s)}\left[\text{CF}_{i}^{(t)}(s,a^{\prime})\right]$
stored in expectation. Because we need to specify x as expectation for the
aggregator function, while expectation is defined upon the returned value of
the aggregate function, we have to put a placeholder of expectation at line 7
before aggregate.

##### Line 17.

Line 17 updates the value of $R^{(t)}_{i}(s,\cdot)$ to
$R^{(t+1)}_{i}(s,\cdot)$ according to Equation 3.1.

##### Line 18.

Line 18 computes $\bm{\pi}_{i}^{(t+2)}(\cdot{\,|\,}s)$ according to Equation
3.1 to replace the old value $\bm{\pi}_{i}^{(t+1)}(\cdot{\,|\,}s)$ stored in
self.strategy.

###  3.4 Loading the Game

LiteEFG is fully compatible with OpenSpiel, _i.e._ LiteEFG supports almost all
games in OpenSpiel444Some games such as Hanabi (Bard et al., 2020) is too big
and the implementation of OpenSpiel does not include infoset, so that LiteEFG
does not support it.. Moreover, for games not implemented by OpenSpiel, users
can write a game description text file alternatively and load it using
LiteEFG.FileEnv. An example of the game file is illustrated in Figure 3, and
the full example can be found in LiteEFG/game_instances/kuhn.game. At the
beginning, the game file displays the parameters of the game, where the
parameter num_players is necessary and other parameters are optional. The next
several lines will include the node information, with an identifier node at
the beginning of the line, and the node’s name goes after it. For different
types of nodes, the additional information should obey the following rules,

  * •

Chance Node ($p(h)=c$): Keyword chance should be placed at first and the
chance events leading by actions will be displayed afterward. The format is
event_name=probability. For instance, in Figure 3, the chance event JQ means
the private card dealt to player 1 and 2 is Jack and Queen individually.

  * •

Player Node ($p(h)\in[N]$): Keyword player should be placed first and $p(h)$
should go after it. Then, the valid actions leading by actions should be
placed afterward. In Figure 3, k, c, f, b is check, call, fold, and bet
individually.

  * •

Leaf (Terminal) Node ($h\in\mathcal{Z}$): Keyword leaf should be placed first
and the utility leading by leading by payoffs go afterward. The format of
utility is $i=\mathcal{U}_{i}(h)$, where $i\in[N]$ is the player index.

After that, the infosets will be described. Each line of infoset description
will be led by the identifier infoset and the name of the infoset will be
placed after it. Then, the nodes in that infoset will be introduced. After the
keyword nodes, the name of those nodes in the infoset will be displayed.

Figure 3: An example of Kuhn Poker (Kuhn, 1950) represented in the game file
format supported by LiteEFG.

###  3.5 Training

Figure 4: Training and evaluation of CFR algorithm defined in Figure 2.

In line 19 and line 20 of Figure 4, we load the environment leduc_poker from
OpenSpiel to LiteEFG. The traverse_type specifies how to traverse the game
tree in each iteration. Currently, LiteEFG supports enumerating all nodes,
external sampling555For each player $i\in[N]$, we traverse the game tree once
according to the following rule. For a node with $p(h)=i$, we visit children
under $(h,a)$ for each $a\in\mathcal{A}_{h}$. For a node with $p(h)\not=i$, we
sample a children $a\sim\pi_{p(h)}(\cdot{\,|\,}s(h))$ and only visit the
children under $(h,a)$. (Lanctot et al., 2009), and outcome
sampling666Whenever meets a node, sample a children
$a\sim\pi_{p(h)}(\cdot{\,|\,}s(h))$ and only visit the children under $(h,a)$.
(Lanctot et al., 2009).

In line 21 and line 22 of Figure 4, we load the computation graph of CFR
defined in Figure 2 into the environment. In line 25, we update the dynamic
part of the graph, so that the update-rule in Equation 3.1 is executed.

CFR algorithm only guarantees the average sequence-form strategy will
converges to the Nash equilibrium in a two-player zero-sum game, that is

$\displaystyle\left(\frac{1}{T}\sum_{t=1}^{T}\mu_{1}^{\pi_{1}^{(t)}},\frac{1}{T}\sum_{t=1}^{T}\mu_{2}^{\pi_{2}^{(t)}}\right)\text{
converging to the Nash equilibrium.}$ (3.2)

Therefore, we need to maintain the average-iterate sequence-form strategy. To
simplify the implementation, LiteEFG provides the function
LiteEFG.Environment.update_strategy, which will automatically maintain the
average-iterate sequence-form strategy (also supports other types of sequence-
form strategy, and the details can be found in the Github repository).

###  3.6 Evaluation

In line 28 of Figure 4, we print the exploitability of the strategy stored at
graph node graph.strategy. "avg-iterate" indicates we want to measure that of
the average-iterate strategy. LiteEFG.Environment.exploitability returns a
vector $\bm{v}$, where each element $v_{i}\geq 0$ indicates how much player
$i$ can improve her utility by deviating to other strategies while the
strategies of other players remain fixed. The sum $\sum_{i\in[N]}v_{i}$ is the
exploitability, which measures the distance to the Nash equilibrium. When the
exploitability is zero, it implies that the current strategies of all players
form a Nash equilibrium.

###  3.7 Debug

To debug the strategy, users can call LiteEFG.OpenSpielEnv.get_strategy to get
a pandas.DataFrame as shown in Figure 5 to display the strategies.

Figure 5: A snapshot of the strategy of player 1 generated by CFR algorithm in
leduc_poker(suit_isomorphism=True) of OpenSpiel. The first column displays the
name of the infoset, and the second to the fourth column is the probability of
choosing fold / call / raise in that infoset. The index of fold / call / raise
in the representation of infoset is 0, 1, 2.

Moreover, the user can also interact with the strategy computed by the
algorithm, by calling LiteEFG.OpenSpielEnv.interact. The interaction is
displayed in Figure 6.

Figure 6: The interaction with strategy generated by CFR in
LiteEFG.OpenSpielEnv.get_strategy.

###  3.8 Various Baselines

Though EFG is developing fast in recent years (Lee et al., 2021; Liu et al.,
2023; Sokota et al., 2023), it is hard to make a fair comparison between
different algorithms. The difficulty is mainly two-fold. Firstly, some
algorithms are not open-sourced and algorithms in EFGs are sensitive to the
hyper-parameters. Therefore, when re-implementing the baseline algorithms,
researchers may not be able to achieve the same performance as the original
paper with the same set of parameters. Moreover, sometimes researchers do not
use OpenSpiel or even not using Python to implement their algorithms.
Secondly, even though some algorithms are open-sourced (Sokota et al., 2023),
they are highly inefficient since they are purely based on Python. As a
result, it takes a lot of computation resources and time to make a fair
comparison with previous results. Therefore, comprehensive and efficient
baselines are important for the whole community, and LiteEFG provides a
variety of baselines. We list the baseline algorithms in Table 1.

Algorithms |  Traverse Type | Reference
---|---|---
Counterfactual Regret Minimization (CFR) |  Full Information | Zinkevich et al. (2007)
Counterfactual Regret Minimization+ (CFR+) |  Full Information | Tammelin et al. (2015)
Discounted Counterfactual Regret Minimization (DCFR) |  Full Information | Brown and Sandholm (2019)
Predictive Counterfactual Regret Minimization (PCFR) |  Full Information | Farina et al. (2021)
External-Sampling Counterfactual Regret Minimization (ES-CFR) |  External Sampling | Lanctot et al. (2009)
Outcome-Sampling Counterfactual Regret Minimization (OS-CFR) |  Outcome Sampling | Lanctot et al. (2009)
Dilated Optimistic Mirror Descent (DOMD) |  Full Information | Lee et al. (2021)
Regularized Dilated Optimistic Mirror Descent (Reg-DOMD) |  Full Information | Liu et al. (2023)
Regularized Counterfactual Regret Minimization (Reg-CFR) |  Full Information | Liu et al. (2023)
Magnetic Mirror Descent (MMD) |  Full Information /
Outcome Sampling | Sokota et al. (2023)
Q-Function Based Regret Minimization (QFR) |  Full Information /
Outcome Sampling | Liu et al. (2024)
Clairvoyant Mirror Descent (CMD) |  Full Information | Wibisono et al. (2022)

Table 1: The baseline algorithms implemented by LiteEFG.

##  4 Benchmark

In this section, we will compare the performance of CFR over LiteEFG and
OpenSpiel. All experiments are computed on Intel(R) Xeon(R) Platinum 8260 CPU
@ 2.40GHz. We compare the performance of the baseline algorithm CFR (Zinkevich
et al., 2007) and CFR+ (Tammelin et al., 2015) in four classical benchmark
games, Liar’s Dice, Leduc Poker (Southey et al., 2005), Kuhn Poker (Kuhn,
1950), and Dark Hex. The results are shown in Figure 7. To make a fair
comparison, we directly call the official C++ implementation of CFR / CFR+ in
OpenSpiel by open_spiel.python.algorithms.cfr._CFRSolver.

We can see that in Figure 7, LiteEFG provides over $100\times$ acceleration in
both games, compared to OpenSpiel. In OpenSpiel, running CFR+ for 100,000
iterations takes about 8 hours in Leduc Poker, while it takes less than 10
minutes with LiteEFG. In the regime of learning in EFGs, the algorithms
usually require extensive hyper-parameter search to enhance the performance
(Lee et al., 2021; Liu et al., 2023; Sokota et al., 2023), which further
enlarges the gap between LiteEFG and OpenSpiel.

CFR (Zinkevich et al., 2007) CFR+ (Tammelin et al., 2015) Figure 7: The
running time of CFR and CFR+ in 4 benchmark games, Liar’s Dice, Leduc Poker,
Kuhn Poker, and Dark Hex. Each algorithm is executed for 100 times in each
game, and for each execution, the algorithm lasts for 100,000 iterations. From
the figure, LiteEFG is about 100$\times$ faster than OpenSpiel.

##  5 Conclusion

In this paper, we propose the new computation framework LiteEFG for solving
EFGs. It is more computationally efficient compared to previous work.
Moreover, users can avoid handling the complex imperfect-information structure
of EFGs by using LiteEFG, since the library will automatically process the
data flow between decision nodes with imperfect information and real nodes in
the game tree. Therefore, LiteEFG would benefit researchers by improving the
efficiency of both implementation of the code and computation for algorithms
in EFGs.

## References

  * Abadi et al. (2016) Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  $\\{$TensorFlow$\\}$: a system for $\\{$Large-Scale$\\}$ machine learning.  In _12th USENIX symposium on operating systems design and implementation (OSDI 16)_ , pages 265–283, 2016.
  * Bard et al. (2020) Nolan Bard, Jakob N Foerster, Sarath Chandar, Neil Burch, Marc Lanctot, H Francis Song, Emilio Parisotto, Vincent Dumoulin, Subhodeep Moitra, Edward Hughes, et al.  The hanabi challenge: A new frontier for ai research.  _Artificial Intelligence_ , 280:103216, 2020.
  * Berner et al. (2019) Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemysław Dębiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, et al.  Dota 2 with large scale deep reinforcement learning.  _arXiv preprint arXiv:1912.06680_ , 2019.
  * Brown and Sandholm (2019) Noam Brown and Tuomas Sandholm.  Solving imperfect-information games via discounted regret minimization.  In _AAAI Conference on Artificial Intelligence (AAAI)_ , 2019.
  * Farina et al. (2021) Gabriele Farina, Christian Kroer, and Tuomas Sandholm.  Faster game solving via predictive blackwell approachability: Connecting regret matching and mirror descent.  In _AAAI Conference on Artificial Intelligence (AAAI)_ , 2021.
  * Kuhn (1950) Harold W Kuhn.  A simplified two-person poker.  _Contributions to the Theory of Games_ , 1(417):97–103, 1950.
  * Lanctot et al. (2009) Marc Lanctot, Kevin Waugh, Martin Zinkevich, and Michael Bowling.  Monte carlo sampling for regret minimization in extensive games.  _Neural Information Processing Systems (NeurIPS)_ , 2009.
  * Lanctot et al. (2019) Marc Lanctot, Edward Lockhart, Jean-Baptiste Lespiau, Vinicius Zambaldi, Satyaki Upadhyay, Julien Pérolat, Sriram Srinivasan, Finbarr Timbers, Karl Tuyls, Shayegan Omidshafiei, Daniel Hennes, Dustin Morrill, Paul Muller, Timo Ewalds, Ryan Faulkner, János Kramár, Bart De Vylder, Brennan Saeta, James Bradbury, David Ding, Sebastian Borgeaud, Matthew Lai, Julian Schrittwieser, Thomas Anthony, Edward Hughes, Ivo Danihelka, and Jonah Ryan-Davis.  OpenSpiel: A framework for reinforcement learning in games.  2019\.
  * Lee et al. (2021) Chung-Wei Lee, Christian Kroer, and Haipeng Luo.  Last-iterate convergence in extensive-form games.  In _Neural Information Processing Systems (NeurIPS)_ , 2021.
  * Liu et al. (2023) Mingyang Liu, Asuman E. Ozdaglar, Tiancheng Yu, and Kaiqing Zhang.  The power of regularization in solving extensive-form games.  In _International Conference on Learning Representations (ICLR)_ , 2023.
  * Liu et al. (2024) Mingyang Liu, Gabriele Farina, and Asuman E. Ozdaglar.  A policy-gradient approach to solving imperfect-information extensive-form games with iterate convergence.  2024\.
  * Mnih et al. (2013) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Silver et al. (2016) David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  _nature_ , 529(7587):484–489, 2016.
  * Silver et al. (2017) David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of go without human knowledge.  _nature_ , 550(7676):354–359, 2017.
  * Sokota et al. (2023) Samuel Sokota, Ryan D’Orazio, J. Zico Kolter, Nicolas Loizou, Marc Lanctot, Ioannis Mitliagkas, Noam Brown, and Christian Kroer.  A unified approach to reinforcement learning, quantal response equilibria, and two-player zero-sum games.  In _International Conference on Learning Representations (ICLR)_ , 2023.
  * Southey et al. (2005) Finnegan Southey, Michael Bowling, Bryce Larson, Carmelo Piccione, Neil Burch, Darse Billings, and Chris Rayner.  Bayes’ bluff: opponent modelling in poker.  In _Conference on Uncertainty in Artificial Intelligence (UAI)_ , 2005.
  * Steiner et al. (2019) Benoit Steiner, Zachary DeVito, Soumith Chintala, Sam Gross, Adam Paske, Francisco Massa, Adam Lerer, Greg Chanan, Zeming Lin, Edward Yang, et al.  Pytorch: An imperative style, high-performance deep learning library.  2019\.
  * Tammelin et al. (2015) Oskari Tammelin, Neil Burch, Michael Johanson, and Michael Bowling.  Solving Heads-Up Limit Texas Hold’em.  In _International Joint Conference on Artificial Intelligence (IJCAI)_ , 2015.
  * Wibisono et al. (2022) Andre Wibisono, Molei Tao, and Georgios Piliouras.  Alternating mirror descent for constrained min-max games.  _Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Zinkevich et al. (2007) Martin Zinkevich, Michael Johanson, Michael H. Bowling, and Carmelo Piccione.  Regret minimization in games with incomplete information.  In _Neural Information Processing Systems (NeurIPS)_ , 2007.