# Hierarchical autoregressive neural networks for statistical systems

Piotr Białas<EMAIL_ADDRESS>Piotr Korcyl<EMAIL_ADDRESS>Tomasz Stebel<EMAIL_ADDRESS>Institute of Applied Computer Science,
Jagiellonian University, ul. Łojasiewicza 11, 30-348 Kraków Poland  Institute
of Theoretical Physics, Jagiellonian University, ul. Łojasiewicza 11, 30-348
Kraków Poland

###### Abstract

It was recently proposed that neural networks could be used to approximate
many-dimensional probability distributions that appear e.g. in lattice field
theories or statistical mechanics. Subsequently they can be used as
variational approximators to asses extensive properties of statistical
systems, like free energy, and also as neural samplers used in Monte Carlo
simulations. The practical application of this approach is unfortunately
limited by its unfavorable scaling both of the numerical cost required for
training, and the memory requirements with the system size. This is due to the
fact that the original proposition involved a neural network of width which
scaled with the total number of degrees of freedom, e.g. $L^{2}$ in case of a
two dimensional $L\times L$ lattice. In this work we propose a hierarchical
association of physical degrees of freedom, for instance spins, to neurons
which replaces it with the scaling with the linear extent $L$ of the system.
We demonstrate our approach on the two-dimensional Ising model by simulating
lattices of various sizes up to $128\times 128$ spins, with time benchmarks
reaching lattices of size $512\times 512$. We observe that our proposal
improves the quality of neural network training, i.e. the approximated
probability distribution is closer to the target that could be previously
achieved. As a consequence, the variational free energy reaches a value closer
to its theoretical expectation and, if applied in a Markov Chain Monte Carlo
algorithm, the resulting autocorrelation time is smaller. Finally, the
replacement of a single neural network by a hierarchy of smaller networks
considerably reduces the memory requirements.

###### keywords:

Variational Autoregressive Neural Networks , Hierarchical Neural Networks ,
Spin Systems , Ising model , Markov Chain Monte Carlo

††journal: Computer Physics Communications

##  1 Introduction

In Markov Chain Monte Carlo simulations the main obstacle in improving the
statistical precision are the autocorrelations between consecutive samples.
The recent idea [1, 2] to alleviate this problem is based on the introduction
of artificial neural networks which can be trained to approximate the target
probability distribution. In such situation the next sample can be generated
independently by the neural network instead of relying on local updates of the
previous sample. If the probability distribution parametrized by the neural
network is sufficiently close to the target probability distribution most of
the proposals would be accepted and the autocorrelations would be reduced.
This algorithm, called Neural Markov Chain Monte Carlo (NMCMC), was studied in
details in Ref. [3] using the two-dimensional Ising model as the testbed. It
was confirmed that the critical slowing down can be much milder than in the
local Metropolis algorithm [4]. Similar conclusions concerning the NMCMC
autocorrelation times, in systems with continuous degrees of freedom, were
reported earlier in Ref. [5].

The work of Ref. [2] was based on a previous observation that autoregressive
neural networks can be applied to estimate extensive observables in
statistical physics, in particular the free energy, which can be directly
related to the Kulback–Leibler divergence, see [6] for details. More recently
that observation, reformulated in terms of normalising flows in order to
accommodate continuous degrees of freedom, was applied to the popular
elementary particle model of scalar $\phi^{4}$ field theory in two-dimensions
in Ref. [7]. This may be considered as true progress since estimating the free
energy of a system with discrete or continuous degrees of freedom is a
notoriously difficult problem for standard Monte Carlo approaches. Another
application of autoregressive neural networks was reported in Ref. [8].

Staying in the context of systems with discrete degrees of freedom, the main
obstacle in scaling the NMCMC approach or the variational estimation of the
free energy, to system sizes relevant from the point of view of particle or
statistical physics is the numerical cost of the approach which dramatically
increases with the system size. Most of the published results used the maximal
lattice of $L\times L=16\times 16$ degrees of freedom. In this work we
describe a mapping of the physical degrees of freedom to the inputs of a set
of neural networks which considerably reduces the numerical cost. We
demonstrate our approach on the two-dimensional Ising model where we simulate
a system of $128\times 128$ spins as a test of validity of our idea. Compared
to the original approach described in Ref. [6] which is valid for any system
with discrete degrees of freedom, our algorithm is applicable only to systems
with nearest-neighbour interactions and it relies on the well-known
Hammersley-Clifford theorem [9, 10] and Markov field property of the Ising
model.

As suggested by Ref. [6], the prohibitive scaling might be eliminated by
replacing fully-connected layers by convolutional ones in the employed
architecture of neural network. Convolutional layers proved their superiority
in applications where the translational symmetry is expected on the level of
input data. This is, however, not the case in the context of autoregressive
networks discussed here, as neurons represent conditional probabilities of
consecutive physical degrees of freedom and there is no symmetry on the level
of conditional probabilities. In particular, the association of physical spins
to consecutive neurons is arbitrary and can even be random. This implies that
the neighbouring neurons do not necessarily correspond to spins which are
close to each other on the physical lattice. This means that convolutional
layers do not enforce the translational symmetry any more in this case. In
fact they may even break it. Fully connected layers also do not enforce
translational symmetry, but at least in theory can be trained to respect it.
Note that in convolutional architectures the number of layers must scale like
$L$ to fully account for all dependencies in conditional probabilities. It was
shown in Ref. [6] that fully connected architectures are much faster than
convolutional layers when actual run times are compared. Hence, in the
following we do not consider convolutional layers as a viable solution for
this particular class of problems.

As announced, we present our idea using the two-dimensional ferromagnetic
Ising model with $N=L^{2}$ spins which are located on a $L\times L$ periodic
lattice. The system is fully characterized by the Hamiltonian operator which
provides the energy of a given configuration of spins
$\mathbf{s}=\\{s^{i}\\}_{i=1}^{N}$

$H(\mathbf{s})=-J\sum_{\langle i,j\rangle}s^{i}\,s^{j}\,,$ (1)

where the sum runs over all neighbouring pairs of lattice sites. We set the
coupling constant $J=1$. Each configuration $\mathbf{s}$ appears with the
Boltzmann probability

$\displaystyle p(\mathbf{s})=\frac{1}{Z}\exp(-\beta H(\mathbf{s}))\,,$ (2)

where $Z$ is the partition function $Z=\sum_{\mathbf{s}}\exp(-\beta
H(\mathbf{s}))$ and the sum is performed over all $2^{N}$ configurations.

##  2 Conditional probabilities and autoregressive networks

In the approach discussed in earlier works [6, 2] a neural network was
proposed to model the probability $p(\mathbf{s})$. To this end $p(\mathbf{s})$
was rewritten as a product of conditional probabilities of consecutive spins,

$p(\mathbf{s})=p(s^{1})\prod_{i=2}^{N}p(s^{i}|s^{1},\dots,s^{i-1})$ (3)

The neural network is used to model the conditional probabilities
$p(s^{i}|\cdot)$ by approximating it by $q_{\theta}(\mathbf{s})$. Then,

$q_{\theta}(\mathbf{s})=q_{\theta}(s^{1})\prod_{i=2}^{N}q_{\theta}(s^{i}|s^{1},\dots,s^{i-1}),$
(4)

where $q_{\theta}$ are quantities returned by the neural network. In the
simplest realization, the neural network has $N$ input and output neurons and
the value of the $i$-th output neuron is interpreted as
$q_{\theta}(s^{i}=+1|s^{1},s^{2},\dots,s^{i-1})$, where we explicitly marked
the dependence of $q$ on the internal parameters $\theta$ of the neural
network. In order to ensure that the conditional probabilities depend only on
the values of the preceding spins one uses autoregressive neural network,
where the neural connections are only allowed to the neurons associated to the
preceding spins.

Knowing the target probability of the spin configuration, $p(\mathbf{s})$ and
its value approximated by the neural network, $q_{\theta}(\mathbf{s})$, the
training of the neural network is performed by minimizing the inverse
Kullback–Leibler (KL) divergence,

$D_{KL}(q_{\theta}|p)=\sum_{\mathbf{s}}\,q_{\theta}(\mathbf{s})\left(\log
q_{\theta}(\mathbf{s})-\log p(\mathbf{s})\right)=E[\log
q_{\theta}(\mathbf{s})-\log p(\mathbf{s})]_{q_{\theta}(\mathbf{s})}.$ (5)

As we usually do not know the normalization factor $Z$, by inserting the
target probability $p$ _without_ the normalization factor $1/Z$ into the
inverse KL divergence $D_{KL}(q_{\theta}|p)$ we get

$\begin{split}F_{q}=E[\log q_{\theta}(\mathbf{s})-\log p(\mathbf{s})-\log
Z)]_{q_{\theta}(\mathbf{s})}=F+D_{KL}(q_{\theta}|p),\end{split}$ (6)

where $F=-\log Z$ is the true free energy which does not depend on the neural
network weights $\theta$, and $F_{q}$ its variational approximation. Hence,
minimization of $D_{KL}(q_{\theta}|p)$ is equivalent to pushing the value of
$F_{q}$ towards the true value $F$ from above. This approach of using neural
networks as a functional ansatz for variational minimization is known as
variational autoregressive networks (VAN) [6]. During the training $F_{q}$ can
approximated by the arithmetic average over $n$ (called _batch size_ , for
which we use the value 1024) samples generated from the distribution
$q_{\theta}$,

$\begin{split}F_{q}\approx\hat{F}_{q}=\frac{1}{n}\sum_{i=1}^{n}\left(\log
q_{\theta}(\mathbf{s}_{i})+\beta
H(\mathbf{s}_{i})\right),\quad\mathbf{s}_{i}\sim q_{\theta}.\end{split}$ (7)

This quantity is minimized during the network training. Its gradient estimator
is calculated as [6]

$\diff{\hat{F}_{q}}{\theta}=\frac{1}{n}\sum_{i=1}^{n}\diff{\log
q_{\theta}(\mathbf{s}_{i})}{\theta}\left[(\log
q_{\theta}(\mathbf{s}_{i})+\beta
H(\mathbf{s}_{i}))-\frac{1}{n}\sum_{j=1}^{n}(\log
q_{\theta}(\mathbf{s}_{j})+\beta
H(\mathbf{s}_{j}))\right],\quad\mathbf{s}_{i}\sim q_{\theta}.$ (8)

Network weights are updated using Adam optimizer[11].

The generation of samples is done via ancestral sampling. We generate spins
one by one starting from $s^{1}$. Once spins $s^{1},\ldots,s^{i-1}$ are
generated they are used as an input to the network to obtain the conditional
probability $q_{\theta}(s^{i}=+1|s^{1},\ldots,s^{i-1})$ according to which the
spin $s^{i}$ is generated and the process is continued. Clearly this requires
calling the network $L^{2}$ times and is the biggest drawback of this
approach. As we are unable to reduce this number, in this contribution we
limit ourselves to the reduction of the workload required by a single call of
the network.

In order to understand the scaling of the numerical cost of the VAN approach
and compare it to our approach we take as our measure two numbers: the
approximate number of trainable parameters of the autoregressive neural
network which is proportional to the memory utilization and the approximate
number of floating point operations needed to generate one configuration of
spins which characterizes the numerical effort. The actual run times, which
depend on the implementation and on the hardware used, are discussed in
Section 5. In order to simplify the discussion we assume that neural networks
are composed of a single fully-connected layer as this is sufficient for our
aims of comparing the two approaches, i.e. including $m_{l}$ layers
approximately multiplies the number of parameters and of floating point
operations by $m_{l}$. We further assume that this layer has same number of
input and output neurons i.e. it looks like the lower layer in Fig. 2. Adding
another layer with smaller number of output neurons does not change our
conclusions about scaling. In what follows we consider large $L\gg 1$ and we
keep only the leading power of $L$ which is denoted with the ”$\approx$” sign.

The original proposal of Ref. [6] assumed a single fully-connected neural
network for all spins and a mask in order to impose autoregressivity, which
corresponds to $L^{4}+L^{2}+L^{2}\approx L^{4}$ trainable parameters. The
first term counts the weights, the second corresponds to the bias on each
neuron and the third counts the parameters in the `PReLU` activation function
used in that Reference. In principle, one could exploit the triangular
(autoregressive) structure of the weight matrix to reduce the number of
parameters by $\frac{1}{2}$. Additionally, one could as well implement another
optimization such as the chessboard factorization [3] and further reduce the
scaling by $L^{2}\rightarrow\frac{1}{2}L^{2}$, giving in total
$\frac{1}{8}L^{4}+L^{2}$. Nevertheless, the number of trainable parameters
would still increase with a leading scaling power of $L^{4}$. As far as the
number of floating point operations is concerned, in order to generate one
configuration of spins we need to invoke the neural network $L^{2}$ times
(number of spins on the lattice) and each invocation has of the order of
$L^{4}$ floating point operations. Hence, in total the number of operations
needed to generate of one configuration scales as $L^{6}$ in the VAN approach.

In the following we propose to reorganize the neural network in such a way as
to keep the same expressivity, but reduce the number of parameters and hence
considerably reduce the numerical cost of the approach allowing us to simulate
systems of size 64 times larger than the original proposal. With much smaller
neural networks our idea considerably reduces memory requirements.

##  3 Boundary and interior spins

Figure 1: Example of hierarchical partitioning for $L=16$. On the first level,
the red boundary $B^{(0)}$ highlighted on the left panel of the figure is
generated with one neural network $\mathcal{N}^{0}$. At the second level of
hierarchy, one neural network of smaller size $\mathcal{N}_{1}$ is used to
consecutively fix four sets of boundaries shown in blue. The example of
$B^{(1),3}$ is highlighted in the middle panel. The surrounding spins
$\partial B^{(1),3}$ are shown in orange. At the third level of hierarchy, one
neural network of even smaller size $\mathcal{N}_{2}$ is used to consecutively
fix sixteen sets of boundary spins marked in green. The example of
$B^{(2),15}$ is highlighted in green on the right panel. The remaining empty
spins corresponding to $I^{k}\equiv B^{(3),k}$, $k=1,\dots,64$ have all the
neighbours fixed and therefore can be generated from a local Boltzmann
distribution with the heatbath algorithm. Figure 2: Scaled down representation
of the architecture of the smallest neural network used to generate green
sublattice (Figure 1 right). Orange nodes correspond to border spins and green
nodes to interior spins. Orange nodes are connected to every node of the
subsequent layer. Green nodes are connected in an autoregressive way. The
propagation of signals starts at the bottom of the figure and goes upward. The
real network used in this case has nine orange input neurons and five green.

We propose to divide the spins into two species: boundary spins and interior
spins, as demonstrated in Fig. 1 (left) for periodic lattice. For non-periodic
lattice the only difference is the shape of the red boundary spins. In the
presented partitioning of a $16\times 16$ lattice we obtain 60 red boundary
spins and four sublattices with 49 white interior spins each. Each sublattice
has boundaries fixed by a subset of red boundary spins. The main advantage of
such a partitioning is the following fact which can be easily proven and is
known as the Hammersley-Clifford theorem: for each sublattice, once the
boundary made of red spins is fixed, the probabilities of interior spins do
not depend on the values of the spins outside the surrounding red boundary.
This allows us to write $p(\mathbf{s})$ as the product

$p(\mathbf{s})=p(\textrm{B}(\mathbf{s}))p(\textrm{I}(\mathbf{s})|\textrm{B}(\mathbf{s}))=p(\textrm{B}(\mathbf{s}))\prod_{a=1}^{4}p(\textrm{I}^{a}(\mathbf{s})|\textrm{B}^{a}(\mathbf{s})),$
(9)

where the operation B, $\textrm{B}(\mathbf{s})$, extracts the boundary spins
from the configuration $\mathbf{s}$, whereas I, $\textrm{I}(\mathbf{s})$ the
interior spins. We explicitly indicated that the probabilities of interior
spins are conditional probabilities depending on the boundary spins. Each
individual factor can be further factorized as

$p(\textrm{B}(\mathbf{s}))=p(s_{B}^{1})\prod_{i=2}^{N_{B}}p(s_{B}^{i}|s_{B}^{1},s_{B}^{2},\dots,s_{B}^{i-1})$
(10)

and

$p(\textrm{I}^{a}(\mathbf{s})|\textrm{B}^{a}(\mathbf{s}))=\prod_{i=1}^{N_{I}}p(s_{I}^{a,i}|s_{I}^{a,1},s_{I}^{a,2},\dots,s_{I}^{a,i-1};\textrm{B}^{a}(\mathbf{s})).$
(11)

At this point we can use _two_ neural networks: one for the boundary spins and
one for the interior spins. The neural network for the boundary spins is of
the same type as in the original NMCMC algorithm, but considerably smaller,
its width is $4L-4$ and hence the number of its trainable parameters scales as
$L^{2}$. The second type of neural network used to parametrize the
probabilities of the interior spins has $(\frac{1}{2}L-1)^{2}$ neurons
corresponding to interior spins (white) and $2L-4$ neurons corresponding to
boundary spins (red), resulting in total of $\approx\frac{1}{4}L^{2}$ of input
neurons. This network has a fully connected part from the boundary input
neurons (orange) to all output neurons, and an autoregressive part from
interior input neurons to interior output neurons (green) (see Figure 2). Note
that the same neural network is needed for all 4 sublattices. This may be
interpreted as introducing a substitute of translational symmetry into the
approach. However, we implement it on the level of physical degrees of freedom
where this symmetry is defined, and not on the level of conditional
probabilities, as would convolutional layers do. Calculating the total number
of trainable weights in this setup we still get the unfavorable scaling of
$L^{4}$ which originates from the neural network used to generate the interior
spins in the sublattices. However, one can easily improve on this, by applying
the division into boundary and interior spins to each sublattice. This will be
described in the next Section.

##  4 Hierarchical Autoregressive Networks (HAN)

Based on the conclusion from the preceding Section, we describe an algorithm
for spin generation based on iterative division of the lattice into smaller
sublattices. Below we describe the procedure and provide a pseudocode which
illustrates the algorithm.

We start again with the division of spins into boundary and interior spins
(see the left panel of Fig. 1). The first set of boundary spins is marked in
red and we generate them as described in the previous Section. We denote these
spins by $B^{(0)}$ where the upper index corresponds to the zeroth step of
consecutive iterations consisting of division and fixing of spins. Similarly,
the corresponding autoregressive neural network trained to fix these spins is
denoted by $\mathcal{N}^{0}$. The spins in $B^{(0)}$ do not depend
conditionally on any other variables. Now, once those spins are fixed, we take
each sublattice and treat it as a new problem with fixed boundaries. Instead
of directly generating all interior spins, we iterate and divide it further
into boundary and remaining spins. The sets of boundary spins $B^{(1),k}$ on
each of the sublattices are marked in blue in the middle panel of Fig. 1. One
example of a new set of spins at the second level of iteration, $B^{(1),3}$,
is highlighted . We propose to use another autoregressive neural network,
called $\mathcal{N}^{1}(\partial B^{(1),k})$ to generate the blue boundary
spins. The main advantage of our approach is that one network
$\mathcal{N}^{1}$ is used to generate all $B^{(1),k},\,k=1,2,3,4$. The
probabilities of blue spins depend conditionally on the surrounding spins, we
denote the latter by $\partial B^{(1)}$. The middle panel in Fig. 1 shows the
example of $B^{(1),3}$ and highlights in orange the surrounding spins
belonging to $\partial B^{(1),k}$. The geometry of the neural network which
must be used in this case has a fully connected part, an autoregressive part,
with the values of the surrounding spins playing the role of external
parameters which is explicitly marked by the argument of the neural network
$\mathcal{N}^{1}$. We show an example of such neural network in Fig. 2.
Repeating this procedure iteratively, we stop when $\frac{3}{4}L^{2}$ boundary
spins (at different levels of iterations) are generated and the remaining
$\frac{1}{4}L^{2}$ interior spins can be set by the heatbath step111In the
heatbath algorithm we draw a spin from probability distribution fixed by the
values of its neighbours: probability of spin $s^{i}$ to be $+1$ or $-1$ is
equal $q(s^{i}=\pm 1,n(s^{i}))=\left[1+\exp\left(\mp\,2\beta\sum_{\ j\in
n(s^{i})}s^{j}\right)\right]^{-1}$, where the sum is performed over neighbours
$n(s^{i})$. (see the the right panel of Fig. 1). The heatbath algorithm may in
principle be replaced by a further level in the hierarchy of neural networks.
However, since that single spin can be generated with exact probability, we
implement that optimization directly. We summarize these steps with the
pseudocode presented in Algorithm 1.

Algorithm 1:  Iterative generation of all spins in a single configuration.

HAN$\beta$$\triangleright$ generate spin configuration $C$ at given $\beta$

define $B^{(0)}$

$B^{(0)}\leftarrow\mathcal{N}^{0}(\textrm{void})$ $\triangleright$ spins
$B^{(0)}$ are set by $\mathcal{N}^{0}$ without any additional dependence

$C\leftarrow B^{(0)}$ $m=1,\dots,\log_{2}L$ $k=1,\dots,4^{m}$

define $B^{(m),k}$ $\triangleright$ one set of blue sites on Fig.1 (middle)

define $\partial B^{(m),k}$ $\triangleright$ set of surrounding spins of
$B^{(m),k}$

$B^{(m),k}\leftarrow\mathcal{N}^{m}(\partial B^{(m),k})$ $\triangleright$
spins $B^{(m),k}$ are set by $\mathcal{N}^{m}$ depending on $\partial
B^{(m),k}$

$C\leftarrow B^{(m),k}$ $\triangleright$ only single spins remain unset, see
white sites on Fig.1 (right)

$I\leftarrow$ HB$(\beta,C)$ $\triangleright$ fix remaining spins with the
heatbath algorithm

$C\leftarrow I$ $\triangleright$ all spins are fixed

return $C$

In order to bound the numerical cost of our approach, we again count the total
number of trainable parameters needed to describe a system of extent
$L=2^{m}=l^{(0)}$, $m>1$. Again, we only keep the leading powers of $L$,
omitting the parameters hidden in the bias and activation functions. The first
step differs from the recursively defined division into boundary and interior
spins as it needs to account for the periodic boundary conditions of the
physical system. We already discussed that step and obtained
$N_{p}^{(0)}=(4L-4)^{2}\approx 16L^{2}=2^{2m+4}$ parameters. The resulting
sublattices have linear extent of
$l^{(1)}=\frac{1}{2}(l^{(0)}-2)=2^{(m-1)}-1$.

If $m>2$, we can continue this process and divide each sublattice by
generating first a new set of boundary spins (blue spins in Fig. 1 middle)
that divide each sublattice in four sublattices. To this end we use a new
neural network. All blue spins in given sublattice depend conditionally on
$4l^{(1)}$ red boundary spins which surrounds them (orange spins in Fig. 1
middle), as already explained previously. We have $2l^{(1)}-1$ blue spins and
hence the fully connected dense layer will have
$N_{p}^{(1)}\approx(6l^{(1)}-1)^{2}\approx 9\ 2^{2m}$ parameters. At each
subsequent recursive partitioning, we find a sublattice of linear size

$l^{(i)}=\frac{1}{2}\big{(}l^{(i-1)}-1\big{)}=2^{m-i}-1$ (12)

and the number of trainable parameters of network $\mathcal{N}^{i}$ is

$N_{p}^{(i)}\approx(6l^{(i)}-1)^{2}\approx 9\ 2^{2(m-i+1)}.$ (13)

The recursion ends when $l^{(i)}=3$, hence $i=m-2$. Gathering the individual
contributions from each iteration we finally obtain

$N_{p}=\sum_{i=0}^{m-2}N_{p}^{(i)}\approx
2^{2m+4}+9\sum_{i=1}^{m-2}2^{2(m-i+1)}\approx 28L^{2}$ (14)

Therefore, the scaling of parameters number with $L^{4}$ present in the
original approach was completely eliminated and replaced with $L^{2}$.

The number of floating point operations for the neural network at $i$th level
of the hierarchy is a product of $N_{p}^{(i)}$ and number of network
invocations (note that for each spin we need one invocation of network):

$\displaystyle C^{(0)}$ $\displaystyle=(2^{m+2}-4)N_{p}^{(0)}\approx 64\
2^{3m},$ (15) $\displaystyle C^{(i)}$
$\displaystyle=4^{i}(2l^{(i)}-1)N_{p}^{(i)}\approx 72\ 2^{3m-i}.$ (16)

Summing up these numbers we obtain the number of floating point operations for
one layer,

$C=C^{(0)}+\sum_{i=1}^{m-2}C^{(i)}\approx 136L^{3},$ (17)

compared to the $\propto L^{6}$ operations in the VAN approach.

The numerical cost of calculating the probability $q(\mathbf{s})$ of a
configuration which has been already generated is given by,

$C_{q}=N_{p}^{(0)}+\sum_{i=1}^{m-2}4^{i}N_{p}^{(i)}\approx 36L^{2}\log_{2}L,$
(18)

since neural networks are invoked only once for each $B^{(i),k}$.

##  5 Numerical results

Figure 3: Time scaling of the VAN and HAN algorithms: we show time of training
for 100 epochs in both approaches as function of system size $L$. Approximate
scaling as $L^{2}$ for HAN and $L^{6}$ for VAN can be seen for sufficiently
large values of $L$ (lines).

In the numerical experiments we used the an adaptation of the `PyTorch`
implementation of autoregressive layers shared by the authors of Ref. [6]. All
neural networks $\mathcal{N}^{i}$ consist of two densely connected layers with
some connections removed. The architecture was schematically shown in Fig. 2.
The neural network used to model the first boundary spins does not have orange
nodes. First and second layers use `PReLU` and `sigmoid` activation functions
respectively.

As a demonstration of the practical impact of our proposal we compared our
approach with the results obtained using standard VAN from Ref. [6]. For a
given size of the spin lattice we measured the time needed to perform 100
epochs of training of a single neural network (VAN) and the hierarchy of
autoregressive networks (HAN) in the new approach. We plot the results
obtained using a single V100 GPU in Fig. 3. We checked that the time of all
additional operations is negligible compared to the neural network
evaluations. We notice that for small values of $L$, $L\leq 16$, VAN approach
is faster and the time of both methods is comparable for $L=32$. For larger
systems the VAN approach scale approximately as $L^{6}$. As for the HAN
approach we clearly see a $L^{2}$ scaling. We checked that these times are
saturated by the multiple invocations of the smallest neural networks. Hence,
the visible $L^{2}$ scaling can be associated with the overhead of numerous
GPU kernel invocations (which indeed increases as $L^{2}$) and hides the
expected numerical effort scaling. Also for very small neural networks where
the number of parameters is smaller then the number of ALUs on the GPU (5120
for the V100) we expect that time for single network invocation independent on
network size. Decreasing the number of GPU calls should further reduce the run
time of the HAN approach. Still, for $L>32$ already at this stage our proposal
yields much faster configuration generation and training.

What we also expect is that the quality of training is much better than in the
variational approach from Ref. [6]. Our recursive model for approximating the
Boltzmann probability distribution of spin configurations should allow to
obtain a variational estimate of the free energy which is much closer to the
exact value $F$. The latter is obtained using known analytical expressions for
Ising model for given $\beta$ and $L$, see Refs. [12, 2].

Figure 4: Relative deviation of variational free energy $\hat{F}_{q}$ from
exact value as a function of the neural network training time on lattices from
$L=8$ up to $L=128$ obtained with the original and hierarchical algorithm. A
large improvement can be observed for the new algorithm, as the precision of
$\hat{F}_{q}$ is below 0.1%.

In order to illustrate this we show in Fig. 4 the relative difference between
the estimated free energy and exact value, $(\hat{F}_{q}-F)/|F|$. We plot it
for several values of $L$ and $\beta=0.44$, where the expected phase
transition is approximately located. For both approaches the same
$\beta$-annealing algorithm was applied: during training the beta was changed
according to $\beta(epoch)=0.44(1-0.996^{epoch})$. We checked that the use or
not of $\beta$-annealing does not effect the final results for the free
energy. $\hat{F}_{q}$ is estimated from the batch during training and is
plotted together with its statistical uncertainty. In the original VAN (left
panel in Fig. 4) the training of the neural network needed to model the system
of size $L=32$ had difficulties to converge and yielded an approximation of
the free energy 1.5% away of the expected result. On the contrary, the HAN
approach offers faster training and the discrepancy between the approximated
free energy and the expected value below 0.1% even for the system of size
$L=128$ (see right panel in Fig. 4).

The results for $\hat{F}_{q}$ shown in Fig. 4 have been obtained with the
$Z_{2}$ symmetry (present in the Ising model) imposed during training as
explained in Ref. [6]: for a sampled configuration $\mathbf{s}$, it requires
replacement
$q_{\theta}(\mathbf{s})\rightarrow\left[q_{\theta}(\mathbf{s})+q_{\theta}(-\mathbf{s})\right]/2$,
where $-\mathbf{s}$ is obtained from $\mathbf{s}$ by flipping all the spins.
As it was mentioned in Ref. [6] inclusion of $Z_{2}$ significantly improves
training efficiency. The effect of this symmetry was then further investigated
in Ref. [3] by the present Authors, where we have shown that imposing the
$Z_{2}$ leads to improvement of training only in the ordered state, namely for
$\beta>0.44$. For the disordered state the imposition of another type of
symmetries for the Ising model (with periodic boundary conditions) turns out
to be beneficial – translational symmetries:

$\mathcal{S}_{i}^{T_{x}}\,\mathbf{s}_{x,y}=\mathbf{s}_{x-i,y},\quad\mathcal{S}_{i}^{T_{y}}\,\mathbf{s}_{x,y}=\mathbf{s}_{x,y-i},\qquad
0\leq i<L,$ (19)

where $\mathbf{s}_{x,y}$ is a spin in the position $(x,y)$ on the $L\times L$
lattice. In Ref. [3] we have shown that imposing $T_{y}$ symmetry implemented
by the replacement,

$q_{\theta}(\mathbf{s})\rightarrow\left[q_{\theta}(\mathbf{s})+q(\mathcal{S}_{1}^{T_{y}}\mathbf{s})+\dots+q(\mathcal{S}_{L-1}^{T_{y}}\mathbf{s})\right]/L,$
(20)

further improves training efficiency.

Figure 5: Effect of symmetries applied during the training: variational free
energy $\hat{F}_{q}$ as a function of the neural network training time on
lattices size $L=32$ obtained when only $Z_{2}$ symmetry is applied and when
one combines $Z_{2}$ and $T_{y}$ symmetries.

We show the effect of this improvement in Fig. 5 where the $T_{y}$ symmetry
was included during the training. With this enhancement we reached the
precision $5\times 10^{-5}$ (red dots) for the system with $L=32$ compared
with $0.002$ for HAN with $Z_{2}$ (green dots) and $0.015$ given by VAN with
$Z_{2}$ (green dots in left panel of Fig. 4). The numerical cost of symmetry
imposition adds up to the cost of configurations generation. Each additional
calculation of $q_{\theta}(\mathbf{s})$ which cost is given by Eq.(18) is
subleading compared to the cost of configurations generation (17) as long as
the number of such additional invocations does not scale as $L$. Translational
symmetries are an exception here, because the number of possible
transformations is given by $L$ in each of the dimensions. Hence, the
imposition of the $T_{y}$ symmetry as indicated by Eq. (20) introduces an
additional power of $L$ to the cost of the symmetry imposition which would
grow like $\sim L^{3}\log_{2}L$. Adding also $T_{x}$ symmetry would further
change scaling into $\sim L^{4}\log_{2}L$.

As was already mentioned, both VAN and HAN approaches can be applied to not
only to the variational approximation of the free energy but also as a
generator of Markov Chain Monte Carlo propositions. Using neural networks
saved after a training consisting of 40000 epochs we generated a Markov chain
of length $N_{MC}$. As was shown in Fig. 4 this is the number of epochs after
which we observe no further improvement in the training for both algorithms
and all system sizes. The consecutive propositions in the Markov chain

$\mathbf{s}_{1},\ldots,\mathbf{s}_{N_{MC}},$ (21)

were accepted with probability

$\min\left(1,\frac{p(\mathbf{s}_{k+1})q_{\theta}(\mathbf{s}_{k})}{p(\mathbf{s}_{k})q_{\theta}(\mathbf{s}_{k+1})}\right).$
(22)

This procedure guarantees that spin configurations are distributed according
to $p$ distribution when $N_{MC}\rightarrow\infty$. When $q_{\theta}\neq p$
some propositions are rejected, leading to repetition of configurations in
MCMC and nonzero autocorrelation. In order to quantify them we define the
normalized autocorrelation function for any observable $\mathcal{O}$ (with
mean $\bar{\mathcal{O}}$ and variance $\operatorname{var}\mathcal{O}$) as

${\Gamma}_{\mathcal{O}}(t)=\frac{1}{\operatorname{var}\mathcal{O}}\frac{1}{N_{MC}}\sum_{i=1}^{N_{MC}}\left(\mathcal{O}(\mathbf{s}_{i})-\bar{\mathcal{O}}\right)\left(\mathcal{O}(\mathbf{s}_{i+t})-\bar{\mathcal{O}}\right),$
(23)

and the integrated autocorrelation time

$\tau^{int}_{\mathcal{O}}=1+2\sum_{i=1}^{t_{max}}\Gamma_{\mathcal{O}}(t).$
(24)

The sum is performed up to the value $t_{max}$, where
${\Gamma}_{\mathcal{O}}(t)$ becomes negative (due to statistical
fluctuations). In what follows we choose energy $H$ given by Eq. (1) as the
operator $\mathcal{O}$ and generate Markov chain of length $N_{MC}=6\cdot
10^{5}$ configurations. We use them to measure $\tau^{int}$ (note that we
skipped $H$ subscript). We choose $\beta=0.44$, where $\tau^{int}$ is maximal
due to the critical slowing down (see Ref. [3] for discussion of this
phenomenon in the normal VAN approach).

Fig. 6 compares the integrated autocorrelation time for VAN and HAN as a
function of system size $L$. In both approaches we observe a rise of
$\tau^{int}$ with $L$. However, it is much faster for the VAN approach. It can
be seen that HAN for system $128\times 128$ achieves similar autocorrelation
time as the VAN algorithm for $32\times 32$. These data did not include the
effects of the imposition of the translational symmetry. Since this symmetry
considerably improves the quality of the training of the neural network, we
expect that the integrated autocorrelation times may be even smaller. However,
due to the increased numerical cost of that symmetry we did not consider it
here.

##  6 Conclusions

In this work we have described a particular scheme of generating spin
configurations using autoregressive neural networks which has a better scaling
property than the reference approach. We have provided estimates of its
numerical cost which showed that it scales with the system size as $L^{3}$ as
opposed to $L^{6}$ of the VAN approach. The main improvement comes from the
fact that we have replaced a single autoregressive neural network of width
$\sim L^{2}$ by a hierarchy of $\log_{2}L-1$ neural networks of decreasing
width. It follows that the number of neural networks weights is much smaller
in the new scheme, thus leading to a much smaller memory consumption. We
expect that, from the practical perspective, the improvement is even larger
because it usually takes less steps of the training algorithm to provide
optimal results for a neural network with smaller number of weights. Our idea
can improve both, the variational method of estimating the free energy of
statistical systems, and the autocorrelation times in the Markov Chain Monte
Carlo simulations. Eventually, we emphasize that the hierarchical structure of
neural networks allowed us to reintroduce translational invariance on the
level of physical degrees of freedom. Further improvements may benefit from
this fact.

Figure 6: Integrated autocorrelation time of the hierarchical algorithm as a
function of the system size $L$.

##  7 Outlook

The following observations may serve as possible future research directions.

First, in this work we have demonstrated that the improved scaling of our
proposal works well for the two-dimensional Ising model. However, it is
possible to generalize our recursive approach to higher dimensions as well as
to systems with more complex interactions. For example spin glasses models
where the coupling constant $J$ depends on the bound can also be simulated. In
a straightforward manner the random values of $J_{ij}$ can be included in the
structure of neural networks on subsequent levels in the hierarchy as
additional external parameters. In this way the conditional dependence on
$J_{ij}$ can be incorporated allowing the neural network to correctly account
for the functional dependence on $J_{ij}$. One may expect that if the training
of such hierarchy of neural networks is successful, it could be performed once
and used for any set of $J_{ij}$ without the need of re-training, in contrast
to the VAN approach [6] or alternative method suggested in Ref. [13].
Moreover, since the neural networks provide global updates, such approach
might be able to overcome the notorious problem of traditional methods of
getting trapped in one of the local minima of the free energy. Recent results
showing that ergodic sampling in a case of complicated energy landscape in the
context of normalizing flows was reported in Refs. [14, 15]. Hence, we expect
significant improvement of the quality of results with the HAN approach
applied to spin glass type of systems.

Second, the approach proposed in this work offers an improved training
strategy, as the training of neural networks for larger spin systems does not
need to be started from random weights. Indeed, small neural networks at the
bottom of the HAN hierarchy can be trained to high accuracy using much smaller
spin systems, their states can be saved and reused for larger spin systems.
This is due to the fact that the conditional probabilities depend only on the
surrounding spins and are independent of the total size of the simulated
system. So in principle when doubling the system linear size $L$, only the two
largest neural networks $\mathcal{N}^{0}$ and $\mathcal{N}^{1}$ have to be
constructed and trained ”from scratch”. The rest of the hierarchy
$\mathcal{N}^{i}$ can be reused from the smaller system. Those networks could
be used as they are, without any further training, but probably better results
can be achieved with some additional training which we would expect to
converge faster.

One may also expect that smaller batch sizes can be used for training. One may
observe that the smaller neural networks at higher iterations of the division
are trained using multiple sublattices of the original configuration. In fact
at iteration $i$, we have $4^{i}$ sublattices which we invoke in the training
process and hence the variance of the gradient is $2^{i}$ smaller than in the
original approach. Hence, we expect that the quality of training of these
inner neural networks will be very high. In the case when such precision is
not necessary, one may try to reduce the batch size.

Eventually, the natural question which arises, whether the hierarchical
generation presented in this work can be generalized to systems with
continuous degrees of freedom. Our preliminary studies show that this can be
done for the scalar $\phi^{4}$ field theory using Normalizing Flows. However,
using the hierarchy of neural networks with dense layers turns out to be not
as beneficial as we have demonstrated in the case of autoregressive networks.
Further studies are needed for more definite conclusions.

## Acknowledgments

Computer time allocation ’plgtmdlangevin2’ and ’plgnnformontecarlo’ on the
Prometheus supercomputer hosted by AGH Cyfronet in Kraków, Poland was used
through the polish PLGRID consortium. T.S. kindly acknowledges support of the
Polish National Science Center (NCN) Grants No. 2019/32/C/ST2/00202 and
2021/43/D/ST2/03375 and support of the Faculty of Physics, Astronomy and
Applied Computer Science, Jagiellonian University Grant No.
2021-N17/MNS/000062. This research was partially funded by the Priority
Research Area Digiworld under the program Excellence Initiative – Research
University at the Jagiellonian University in Kraków.

## References

  * Albergo et al. [2019] M. S. Albergo, G. Kanwar, P. E. Shanahan,  Flow-based generative models for markov chain monte carlo in lattice field theory,  Phys. Rev. D 100 (2019) 034515. URL: https://link.aps.org/doi/10.1103/PhysRevD.100.034515. doi:10.1103/PhysRevD.100.034515.
  * Nicoli et al. [2020] K. A. Nicoli, S. Nakajima, N. Strodthoff, W. Samek, K.-R. Müller, P. Kessel,  Asymptotically unbiased estimation of physical observables with neural samplers,  Phys. Rev. E 101 (2020) 023304. doi:10.1103/PhysRevE.101.023304. arXiv:1910.13496.
  * Białas et al. [2021] P. Białas, P. Korcyl, T. Stebel,  Analysis of autocorrelation times in Neural Markov Chain Monte Carlo simulations (2021). arXiv:2111.10189.
  * Metropolis et al. [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller,  Equation of state calculations by fast computing machines,  The Journal of Chemical Physics 21 (1953) 1087–1092. URL: https://doi.org/10.1063/1.1699114. doi:10.1063/1.1699114.
  * Kanwar et al. [2020] G. Kanwar, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, S. Racanière, D. J. Rezende, P. E. Shanahan,  Equivariant flow-based sampling for lattice gauge theory,  Phys. Rev. Lett. 125 (2020) 121601. doi:10.1103/PhysRevLett.125.121601. arXiv:2003.06413.
  * Wu et al. [2019] D. Wu, L. Wang, P. Zhang,  Solving Statistical Mechanics Using Variational Autoregressive Networks,  Phys. Rev. Lett. 122 (2019) 080602. doi:10.1103/PhysRevLett.122.080602. arXiv:1809.10606.
  * Nicoli et al. [2021] K. A. Nicoli, C. J. Anders, L. Funcke, T. Hartung, K. Jansen, P. Kessel, S. Nakajima, P. Stornati,  Estimation of thermodynamic observables in lattice field theories with deep generative models,  Phys. Rev. Lett. 126 (2021) 032001. URL: https://link.aps.org/doi/10.1103/PhysRevLett.126.032001. doi:10.1103/PhysRevLett.126.032001.
  * Luo et al. [2021] D. Luo, Z. Chen, K. Hu, Z. Zhao, V. M. Hur, B. K. Clark,  Gauge Invariant Autoregressive Neural Networks for Quantum Lattice Models (2021). arXiv:2101.07243.
  * J. M. Hammersley [1971] P. C. J. M. Hammersley, Markov fields on finite graphs and lattices, 1971. URL: http://www.statslab.cam.ac.uk/~{}grg/books/hammfest/hamm-cliff.pdf.
  * Clifford [1990] P. Clifford,  Markov random fields in statistics,  in: Disorder in Physical Systems. A Volume in Honour of John M. Hammersley, Clarendon Press, 1990\.
  * Kingma and Ba [2014] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, 2014. URL: https://arxiv.org/abs/1412.6980. doi:10.48550/ARXIV.1412.6980.
  * Ferdinand and Fisher [1969] A. E. Ferdinand, M. E. Fisher,  Bounded and inhomogeneous ising models. i. specific-heat anomaly of a finite lattice,  Phys. Rev. 185 (1969) 832–846. URL: https://link.aps.org/doi/10.1103/PhysRev.185.832. doi:10.1103/PhysRev.185.832.
  * McNaughton et al. [2020] B. McNaughton, M. V. Milošević, A. Perali, S. Pilati,  Boosting Monte Carlo simulations of spin glasses using autoregressive neural networks,  Phys. Rev. E 101 (2020) 053312. doi:10.1103/PhysRevE.101.053312. arXiv:2002.04292.
  * Hackett et al. [2021] D. C. Hackett, C.-C. Hsieh, M. S. Albergo, D. Boyda, J.-W. Chen, K.-F. Chen, K. Cranmer, G. Kanwar, P. E. Shanahan,  Flow-based sampling for multimodal distributions in lattice field theory (2021). arXiv:2107.00734.
  * Albergo et al. [2022] M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, G. Kanwar, S. Racanière, D. J. Rezende, F. Romero-López, P. E. Shanahan, J. M. Urban,  Flow-based sampling in the lattice Schwinger model at criticality (2022). arXiv:2202.11712.


# Electrical analogue of one-dimensional and quasi-one-dimensional Aubry-
André-Harper lattices

Sudin Ganguly<EMAIL_ADDRESS>Department of Physics, School of Applied
Sciences, University of Science and Technology Meghalaya, Ri-Bhoi-793101,
India  Santanu K. Maiti<EMAIL_ADDRESS>Physics and Applied
Mathematics Unit, Indian Statistical Institute, 203 Barrackpore Trunk Road,
Kolkata-700108, India

###### Abstract

The present work discusses the possibility to realize correlated disorder in
electrical circuits and studies the localization phenomena in terms of two-
port impedance. The correlated disorder is incorporated using the Aubry-André-
Harper (AAH) model. One-dimensional and quasi-one-dimensional AAH structures
are explored and directly mapped with their tight-binding analogues.
Transitions from the high-conducting phase to the low-conducting one are
observed for the circuits.

###### pacs:

72.80.Vp, 72.25.-b, 73.23.Ad,

The Aubry-André-Harper (AAH) model aah1 ; aah2  represents a classic example
of a one-dimensional (1D) quasi-crystal, possesses several intriguing
features. In the nearest-neighbor tight-binding (TB) framework, the 1D AAH
model with an incommensurate potential exhibits a sharp localization-
delocalization transition, where all the eigenstates are delocalized below a
critical point, while all of them are completely localized beyond that
critical point. Due to the incommensurate potential, this 1D model shows a
gapped and fractal-like energy spectrum. Beyond the minimal nearest-neighbor
TB model, energy-dependent mobility edges have also been predicted
analytically biddle  in the 1D AAH chain. Such localization transition and
mobility edges have also been found in coupled AAH chains rossi . Recently, 1D
AAH quasi-crystal has been realized experimentally using waveguides kraus  by
Kraus et al. The authors in their work showed that the edge states of the
fabricated photonic quasi-crystal are topologically nontrivial. Owing to such
striking properties and several others, this model has been investigated
widely in many contexts over more than three decades ganeshprl110 ;
ganeshprl114 ; ganeshprb91 ; harter ; purka ; longhi ; wang ; yoo ; lin .

Recently, it has been shown that various topological states that are difficult
to observe in condensed matter experiments, can be simulated with electric
circuits thomale1 ; thomale2 ; helbig . Even though electric circuits
represent classical systems, with a proper choice of the reactive elements,
the corresponding admittance matrix becomes equivalent to the tight-binding
Hamiltonian thomale1 ; thomale2 ; helbig ; jalil-com-phys ; broy . By means of
circuits, the energetics and topological phases of various physical systems
have been investigated in recent years, such as Su-Schrieffer-Heeger (SSH)
model thomale2 , Weyl semimetal jalil-com-phys , Chern and quantum spin Hall
insulators broy , topological Anderson insulators zqzhang , breathing kagome
and pyrochlore lattices ezawaprb98 , and many others zhao ; nakata ; imhof ;
haenel ; zhu ; song ; jalilprb . Most of the aforementioned works focused on
the study of topological phases based on the close correspondence between such
electric circuits and TB models. The motivation of this work is twofold –
simulation of a TB AAH system using an electrical circuit and inspection of
localization phenomena.

First, we construct electric circuits comprising inductors and capacitors (see
Fig. 1) that describe TB AAH systems. The cosine modulation of usual AAH TB
site energies is incorporated into the circuit by connecting different values
of capacitors at the nodes of the circuit.

Figure 1: (Color online.) (a) Schematic diagram of a 1D quasiperiodic TB chain
consisting of eight sites. $t_{1}$ and $t_{2}$ are the nearest-neighbor and
next nearest-neighbor hopping integrals, respectively. (b) Schematic circuit
analogous to 1D quasiperiodic TB chain. (c) Blue capacitors
$\left(C_{N}\text{'s}\right)$ connect the nearest-neighbor nodes. (d) Red
capacitors $\left(C_{NN}\text{'s}\right)$ connect the next nearest-neighbor
nodes. (e) The grounded connection of each node $n$.

Primarily, we focus on designing circuits that are analogous to the TB 1D
chains with nearest-neighbor (NN) and next nearest-neighbor (NNN) connections
and then a two-stranded ladder network. We detect any localization behavior
present in those circuits by computing a two-port impedance (TPI) 2port .

The notable features of this Letter are: (i) realization of correlated
disordered systems with electrical components, (ii) direct mapping of NN and
NNN AAH TB chains and two-stranded AAH ladder with electrical circuits, and
(iii) exact correspondence of admittance spectra of electrical circuits with
energy spectra of TB AAH systems. Our analysis can be implemented to any other
such fascinating correlated systems.

The TB Hamiltonian modeled on a 1D chain (Fig. 1(a)) within a non-interacting
electron picture considering both NN and NNN hoppings can be written as rossi
,

$H=\sum\limits_{n}\epsilon_{n}c_{n}^{\dagger}c_{n}+t_{1}\sum\limits_{\langle
nm\rangle}c_{n}^{\dagger}c_{m}+t_{2}\sum\limits_{\langle\langle
nm\rangle\rangle}c_{n}^{\dagger}c_{m}.$ (1)

Here $t_{1}$ is the NN hopping integral and $\langle\rangle$ represents the NN
sites of the 1D chain, while $t_{2}$ is the NNN hopping strength and
$\langle\langle\rangle\rangle$ denotes the NNN sites of the 1D chain.
$\epsilon_{n}$ is the on-site potential at site $n$. The AAH disorder is
introduced through the on-site potential and it is aah2

$\epsilon_{n}=W\cos{(2\pi bn+\phi_{\nu})}$ (2)

where $W$ is the modulation strength, $b$ is an irrational number and it is
chosen as $b=(\sqrt{5}-1)/2$, $n$ is the site index, and $\phi_{\nu}$ is the
AAH phase factor.

To map the TB 1D lattice model (Fig. 1(a)) that satisfies Eq. 1, we design an
electrical circuit which is given in Fig. 1(b). Kirchhoff’s law at node $n$ in
the given LC circuit reads as broy

$\dot{I}_{n}=\sum_{m}C_{nm}\left(\ddot{V}_{n}-\ddot{V}_{m}\right)+C_{n}\ddot{V}_{n}+\frac{1}{L_{n}}V_{n}.$
(3)

Here $\dot{V}=\frac{dV}{dt}$. $I_{n}$ and $V_{n}$ are the current and voltage
at node $n$, respectively. The sum over $m$ is taken for the first and second
nearest-neighbor nodes. $C_{nm}$ is the capacitor connected between nodes $n$
and $m$. In the second term of Eq. 3, the capacitor $C_{n}$ and the inductor
$L_{n}$ (in the third term) are connected between the node $n$ and ground.

Following the Fourier transformation of Eq. 3, the relationship between the
current and voltage at frequency $\omega$ becomes

$I_{n}(\omega)=\sum_{m}J_{nm}(\omega)V_{m}(\omega).$ (4)

Here $J_{nm}$ is known as the admittance matrix and it becomes

$J_{nm}=j\omega\left[C_{nm}+\left(C_{n}+\left(\sum_{m}C_{nm}-\frac{1}{\omega^{2}L_{n}}\right)\right)\delta_{nm}\right].$
(5)

A one-to-one correspondence can be established between the admittance matrix
$J_{nm}$ and the TB Hamiltonian (Eq. 1) with

$\displaystyle t_{1(2)}$ $\displaystyle=$ $\displaystyle j\omega C_{nm}$ (6)
$\displaystyle\epsilon_{n}$ $\displaystyle=$ $\displaystyle
j\omega\left(C_{n}+\left(\sum_{m}C_{nm}-\frac{1}{\omega^{2}L_{n}}\right)\right).$
(7)

Apart from the term $j\omega$, the NN hopping integral $t_{1}$ can be
identified with the capacitor $C_{N}$ (Fig. 1(b)) and the NNN hopping integral
$t_{2}$ with $C_{NN}$ (Fig. 1(c)). Setting the frequency of the input as
$\omega=1/\sqrt{L_{n}\sum_{m}C_{nm}}$, the term
$\left(\sum_{m}C_{nm}-\frac{1}{\omega^{2}L_{n}}\right)$ in Eq. 7 becomes zero.
Then the capacitor $C_{n}$ (connected between the $n$-th node and the ground)
can be used to incorporate the AAH disorder into the admittance matrix.

To study the localization behavior of the circuit, the simplest experimentally
measurable quantity is the two-port impedance $Z_{nm}$ between nodes $n$ and
$m$. $Z_{nm}$ is defined as 2port

$Z_{nm}=\frac{V_{n}-V_{m}}{I}$ (8)

where $V_{n}-V_{m}$ is the voltage difference between the nodes $n$ and $m$.
$I$ is the magnitude of current $I=I_{n}=-I_{m}$, that is the current $I$
flows into node-$n$ and leaves node-$m$.

In order to determine $Z_{nm}$, we need to express the potentials in terms of
the input current $I$ and for that, the admittance matrix (Eq. 5) has to be
inverted. To do so, first, we write the spectral form of $J_{nm}$ as

$J_{nm}=\sum_{p}j_{p}\psi^{*}_{p,n}\psi_{p,m}$ (9)

where $j_{p}$ is the $p$-th eigenvalue of the admittance matrix and
$\psi_{p,n}$ is the $p$-th eigenfunction at node $n$. With this, the
regularized inverse of the admittance matrix, known as the circuit Green’s
function, can be written as

$G_{nm}=\sum_{p,j_{p}\neq 0}j^{-1}_{p}\psi^{*}_{p,n}\psi_{p,m}.$ (10)

The two-port impedance then simplifies to

$\displaystyle Z_{nm}$ $\displaystyle=$ $\displaystyle
G_{nn}+G_{mm}-G_{nm}-G_{mn}$ (11) $\displaystyle=$
$\displaystyle\sum_{p}j^{-1}_{p}\lvert\psi_{p,n}-\psi_{p,m}\rvert^{2}.$

First, we consider the circuit which is analogous to a 1D TB NN AAH chain. For
this case, we do not consider any connection with the capacitors $C_{NN}$’s.
Before we discuss the results, let us specify the values of the capacitors and
inductors. We consider the number of nodes in the circuit as $N=100$ and set
$C_{N}=1\,\mu$F. The inductor $L_{n}$ that is connected between $n$-th node
and ground is set at $L_{n}=L=1\,$mH except at the extreme two nodes, namely
nodes 1 and 100. The inductors between the ground and these two nodes are
fixed at $L_{1}=L_{100}=2L$. The frequency of the input is fixed at
$\omega=1/\sqrt{2LC_{N}}\,(i\neq 1,100)$. The capacitor $C_{n}$ (Fig. 1(e)) is
considered in the $\mu$F range and the magnitude is chosen according to Eq. 2.
For instance, the grounded capacitors for a disorder strength $W=1$ (with
$\phi_{\nu}=0$) assume the values, $C_{1}=-0.7374\,\mu$F,
$C_{2}=0.0874\,\mu$F, $C_{3}=0.6084\,\mu$F, and so on. In order to observe the
localization behavior in the circuit, the two-port impedance is the measurable
quantity. Now, from the mathematical expression of two-port impedance, one
sees that the term in the denominator is the eigenvalues of the admittance
matrix. If any one of the eigenvalues becomes close to zero, then the two-port
impedance is ought to be very large, irrespective of the nature of the
corresponding eigenstate. Thus, to extract the actual nature of a state, that
is whether the state is localized or not, we attach another capacitor
$C_{\text{offset}}$ (not shown in Fig. 1) between each node and the ground, in
parallel to $C_{n}$.

Figure 2: (Color online.) Admittance spectrum as a function of Aubry phase
$\phi_{\nu}$ considering NN connection. Number of nodes $N=100$. The AAH
disorder strength is fixed at $W=0.5$. Other circuit parameters are described
in the text. The nodes within the red, orange, and black circles are picked up
to compute the two-port impedance in Fig. 3.

This simply shifts all the eigenvalues well above the zero line and hence
makes it possible to get the true localized behavior of the eigenstates. We
set $C_{\text{offset}}=10\,$mF. The effect of such a connection can readily be
observed in the admittance spectra as shown in Fig. 2.

In Fig. 2 we plot the eigenvalues $j_{n}$ of the admittance matrix as a
function of AAH phase $\phi_{\nu}$. As mentioned earlier, the structures of
the admittance and the TB Hamiltonian matrices are identical, apart from the
term $j\omega$ (Eq. 5). The eigenvalues $j_{n}$ of the admittance matrix are
expressed in units of $\Omega^{-1}$ and $\phi_{\nu}$ in units of $\pi$. The
modulation strength is fixed at $W=0.5$. As mentioned above,
$C_{\text{offset}}$ shifts all the eigenvalues towards a positive value and
none of the eigenvalues is close to zero. The spectrum is divided into three
branches and they are almost constant with $\phi_{\nu}$. A few modes are also
seen to cross from one branch to another through the gaps with $\phi_{\nu}$.
Overall, the behavior of the admittance spectrum as a function of $\phi_{\nu}$
is identical to the energy spectrum as a function of $\phi_{\nu}$, computed
for 1D tight-binding Aubry chain by Kraus et al. kraus , and thus, we can
claim that our circuit setup is correct.

Now, let us look at the localization behavior of the circuit. Here we consider
two particular modes in the gap, namely the modes at $\phi_{\nu}=0.48\pi$ (the
red hollow circle in Fig. 2) and $\phi_{\nu}=0.68\pi$ (the orange hollow
circle in Fig. 2), having the same eigenvalues. We compute the two-port
impedance for these two modes as a function of node index as shown in Fig.
3(a).

Figure 3: (Color online.) Visualizing edge states. Two port impedance $|Z|$ as
a function of node index. One port is fixed at node 50. (a) The red curve is
for $\phi_{\nu}=0.48\pi$ and the considered eigenmode is chosen from Fig 2,
marked with red circle. The orange curve is for $\phi_{\nu}=0.68\pi$ and the
considered eigenmode is chosen from Fig 2, marked with orange circle. Both
modes have the same eigenvalue. (b) The black curve is for
$\phi_{\nu}=0.48\pi$ and the considered eigenmode is chosen from Fig 2, marked
with black circle.

We fixed one port at node 50 and the other one is taken through all the nodes
in order to compute the two-port impedance. The results for
$\phi_{\nu}=0.48\pi$ and $\phi_{\nu}=0.68\pi$ are displayed in red and orange
colors, respectively. In Fig. 3(a), we see that the two-port impedance $\lvert
Z\rvert$ is maximum at the extreme left node, namely node 1 for
$\phi_{\nu}=0.48\pi$ and abruptly decreases to zero from the left side to the
right of the 1D electric circuit. On the other hand, for $\phi_{\nu}=0.68\pi$,
we observe a complete mirror-symmetric feature of the previous case. The
maximum impedance is now observed at node 100, which is at the extreme right
of the 1D circuit and then gradually decreases from the right side to the
left. We also choose a mode from the bulk band as shown by the black circle in
Fig. 2 and computed the two-port impedance. The corresponding result is shown
in Fig. 3(b). In the given case, $\lvert Z\rvert$ is about an order of
magnitude less than the previous two cases, indicating that all the nodes are
well extended for the chosen bulk mode. Overall, the behavior of two-port
impedance with node index is very much consistent with the established
localization behavior of 1D TB AAH chain kraus  as a function of $\phi_{\nu}$.

Next, we vary the disorder strength via the grounded capacitors $C_{n}$s and
study the behavior of two-port impedance. Here we set $\phi_{\nu}=0$ and all
other circuit parameters are taken same as mentioned earlier. In the density
plot of Fig. 4, we show the behavior of the natural log of two-port impedance
as functions of disorder strength and eigenvalues of admittance matrix
(measured in $\Omega^{-1}$). Here the two-port impedance

Figure 4: (Color online.) Density plot of the natural log of two-port
impedance as functions of modulation strength $W$ and eigenvalues of
admittance matrix $j_{n}$, considering NN connection. Number of nodes $N=100$.
The colorbar denotes the values of $\ln{\lvert Z\rvert}$. The green vertical
line at $W=2C_{N}$ denotes the sharp transition from the highly conducting
region to the low conducting one. Thus the critical $W$ is $W_{c}=2C_{N}$.

at each of the eigenmodes is computed by keeping one port at node-1 and
varying the other port at all the other nodes and then taking the maximum of
$\lvert Z\rvert$. The values in the colorbar denote the natural log of two-
port impedance, where lower values are denoted with dark color and higher
values with bright ones. The eigenspectrum of the admittance matrix is divided
into three branches as expected. Below $W=2C_{N}$, the computed two-port
impedance for all the modes is vanishingly small, as is clearly seen in Fig.
4. Beyond $W=2C_{N}$, the eigenvalue spectrum becomes brighter, indicating
that the TPI for all the modes is much higher than that in $W<2C_{N}$ region.
Therefore, a sharp transition occurs at $W=2C_{N}$ (critical $W$ is
$W_{c}=2C_{N}$) from a highly conducting zone (vanishingly small $\lvert
Z\rvert$) to a low conducting one (relatively large $\lvert Z\rvert$). Such a
sharp transition in the context of localization has already been studied in 1D
TB AAH chain ganeshprl114 .

Now, we bring in the capacitor $C_{NN}$ to make the circuit analogous to a 1D
NNN TB AAH chain. We fix the capacitor $C_{NN}=0.25\,\mu$F. The inductor
$L_{n}$ that is connected between $n$-th node and ground is set at
$L_{n}=L=1\,$mH except at the extreme four nodes, namely nodes 1, 2, 99, and
100. We fix $L_{1}=L_{100}=2L$ and
$L_{2}=L_{99}=2\left(C_{N}+C_{NN}\right)L/\left(2C_{N}+C_{NN}\right)$. The
frequency of the input is fixed at
$\omega=1/\sqrt{2L\left(C_{N}+C_{NN}\right)}$. The rest of the circuit
parameters are the same as mentioned earlier. The density plot of the TPI,
shown in Fig. 5, is computed following the prescription described above in
Fig. 4. The spectral nature of the eigenvalues of the admittance matrix with
disorder strength is identical to the behavior of the eigenvalue spectrum as a

Figure 5: (Color online.) Density plot of the natural log of two-port
impedance as functions of modulation strength $W$ and eigenvalues of
admittance matrix $j_{n}$, considering NNN connection. Number of nodes
$N=100$. The colorbar denotes the values of $\ln{\lvert Z\rvert}$.

function of disorder strength for the 1D NNN TB AAH chain. Here also we
observe a transition from a conducting region to an insulating one. However,
there is no such sharp transition as in the case of NN case (Fig. 4). Rather,
the transition is admittance dependent. It is important to note that in the 1D
AAH tight-binding chain with higher order hopping terms, there exists energy
dependent mobility edge, which separates the localized wave functions from the
delocalized ones. We also have a similar situation in the present case –
admittance dependent mobility edge, which separates the highly conducting
region from the low conducting one.

Finally, we design a two-stranded ladder network electrically. The schematic
diagram for the circuit is shown in Fig. 6(a). The upper and lower strands are
coupled vertically through the capacitor $C_{V}$ (Fig. 6(b)).

Figure 6: (Color online.) (a) Schematic diagram of two-stranded ladder network
having $N$ nodes in each strand. (b) Blue capacitors $C_{V}$’s connect the two
strands vertically. (c) Red capacitors $C_{H}$’s connect the nearest-neighbor
nodes horizontally in each strand. (d) Cyan capacitors $C_{D}$’s connect the
nodes diagonally. (e) The grounded connection of each node $n$.

In both the strands, the neighboring nodes (marked with green circles) are
connected through the capacitor $C_{H}$, denoted with red color (Fig. 6(c)).
The crossed nodes (viz, nodes 1 and 5, nodes 2 and 4, etc.) are connected
through the capacitor $C_{D}$ (Fig. 6(d)). The $n$-th node is connected to the
ground via a parallel $LC$ circuit, as shown in Fig. 6(e).

The chosen circuit parameters are as follows. The number of nodes fixed at
each strand is $N=100$. The frequency of the input is fixed at
$\omega=1/\sqrt{\left(2C_{H}+2C_{D}+C_{V}\right)L}$. The inductor $L_{n}$ that
is connected between $n$-th node and ground is set at $L_{n}=L=1\,$mH ($n\neq
1,100,101,200$). The inductors between the ground and these four nodes are
fixed at
$L_{n}=\left(2C_{H}+2C_{D}+C_{V}\right)L/\left(C_{H}+C_{D}+C_{V}\right)\,(n=1,100,101,200)$.
The capacitor $C_{n}=C_{n+100}$ (Fig. 1(e)) is considered in the $\mu$F range
and the magnitude is chosen according to Eq. 2 as before. In addition to that,
we attach another capacitor $C_{\text{offset}}=10\,$mF (not shown in Fig. 6)
between each node and ground, in parallel to $C_{n}$ to introduce a shift in
the admittance spectrum well above the zero line due to the fact mentioned
earlier.

With all the said circuit parameters, we show the behavior of the natural log
of two-port impedance as functions of disorder strength and eigenvalues of the
admittance matrix in Fig. 7. The computed two-port impedance for a particular
eigenmode is chosen by considering all the possible two ports in the present
circuit setup and then we take the maximum impedance among them. Here, the
natural log of the maximum impedance is plotted for the density plot. The
color convention for the colorbar is same as before.

Figure 7: (Color online.) Density plot of the natural log of two-port
impedance as functions of modulation strength $W$ and eigenvalues of
admittance matrix $j_{n}$ for two-stranded ladder network. Number of nodes at
each strand $N=100$. The colorbar denotes the values of $\ln{\lvert Z\rvert}$.
The magenta and green vertical lines denote are drawn at two critical points
as described in the texts.

In the present case, we identify two critical points. The first one is
$W_{c1}=2\left(C_{H}-C_{D}\right)=1.6\,\mu$F (shown by the dark magenta
vertical line), below which all the modes are highly conducting. The second
one is $W_{c2}=2\left(C_{H}+C_{D}\right)=2.4\,\mu$F, beyond which all the
modes are poorly conducting in nature (shown by the dark green vertical line).
Within the range $W_{c1}<W<W_{c2}$, there is a mixed phase zone, where the
highly conducting modes and poorly conducting modes coexist. Such a feature is
also in good agreement with the localization behavior of the two-stranded
ladder network in the TB framework rossi .

To conclude, we have proposed a way to realize AAH disorder in electrical
circuits. One-dimensional and two-stranded ladder networks have been
considered for the purpose. We have shown that the AAH disorder strength and
the phase can be controlled by tuning the reactive elements of the circuits.
Like the inverse participation ratio (IPR), which is one of the measures of
the localization phenomena in TB systems, the two-port impedance of the
electrical circuits considered in the present work can serve the same purpose.
Specifically, we have shown that for 1D NN circuit, the behavior of two-port
impedance exhibits a sharp transition from a highly conducting region to a
poorly conducting region. We have also observed admittance dependent mobility
edge in 1D NNN circuit, which separates the high-conducting region from the
low-conducting one. Finally, for the two-stranded ladder network, we have
found two critical points, below one of the critical points, all the modes are
highly conducting, and beyond the other critical point, all of them are poorly
conducting. In between these two critical points, both the low and high
conducting modes coexist. All the observations have been carried out based on
two-port impedance. We strongly believe that the present analysis provides a
direct mapping of AAH lattices with electrical circuits.

## References

  * (1) P. G. Harper, Single band motion of conduction electrons in a uniform magnetic field, Proc. Phys. Soc. London A 68, 874 (1955).
  * (2) S. Aubry and G. André, Analyticity breaking and Anderson localization in incommensurate lattices, Ann. Israel Phys. Soc. 3, 133 (1980).
  * (3) J. Biddle and S. Das Sarma, Predicted Mobility Edges in One-Dimensional Incommensurate Optical Lattices: An Exactly Solvable Model of Anderson Localization, Phys. Rev. Lett. 104, 070601 (2010).
  * (4) M. Rossignolo and L. Dell’Anna, Localization transitions and mobility edges in coupled Aubry-André chains, Phys. Rev. B 99, 054211 (2019).
  * (5) Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin, and O. Zilberberg, Topological States and Adiabatic Pumping in Quasicrystals, Phys. Rev. Lett. 109, 106402 (2012).
  * (6) S. Ganeshan, Kai Sun, and S. Das Sarma, Topological Zero-Energy Modes in Gapless Commensurate Aubry-André-Harper Models, Phys. Rev. Lett. 110, 180403 (2013).
  * (7) S. Ganeshan, J. H. Pixley, and S. Das Sarma, Nearest Neighbor Tight Binding Models with an Exact Mobility Edge in One Dimension, Phys. Rev. Lett. 114, 146601 (2015).
  * (8) S. Ganeshan and S. Das Sarma, Constructing a Weyl semimetal by stacking one-dimensional topological phases, Phys. Rev. B 91, 125438 (2015).
  * (9) A. K. Harter, T. E. Lee, and Y. N. Joglekar, ${\mathcal{PT}}$-breaking threshold in spatially asymmetric Aubry-André and Harper models: Hidden symmetry and topological states, Phys. Rev. A 93, 062101 (2016).
  * (10) A. Purkayastha, S. Sanyal, A. Dhar, and M. Kulkarni, Anomalous transport in the Aubry-André-Harper model in isolated and open systems, Phys. Rev. B 97, 174206 (2018).
  * (11) S. Longhi, Topological Phase Transition in non-Hermitian Quasicrystals, Phys. Rev. Lett. 122, 237601 (2019).
  * (12) Y. Yoo, J. Lee, and B. Swingle, Nonequilibrium steady state phases of the interacting Aubry-André-Harper model, Phys. Rev. B 102, 195142 (2020).
  * (13) Y. Wang, Xu Xia, L. Zhang, H. Yao, S. Chen, J. You, Q. Zhou, and X.-J. Liu, One-Dimensional Quasiperiodic Mosaic Lattice with Exact Mobility Edges, Phys. Rev. Lett. 125, 196604 (2020).
  * (14) Q. Lin, T. Li, L. Xiao, K. Wang, W. Yi, and P. Xue, Topological Phase Transitions and Mobility Edges in Non-Hermitian Quasicrystals, Phys. Rev. Lett. 129, 113601 (2022).
  * (15) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Topolectrical circuit realization of topological corner modes, Nat. Phys. 14, 925 (2018).
  * (16) C. H. Lee, S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, and R. Thomale, Topolectrical Circuits, Commun. Phys. 1, 39 (2018).
  * (17) T. Helbig, T. Hofmann, C. H. Lee, R. Thomale, S. Imhof, L. W. Molenkamp, and T. Kiessling, Band structure engineering and reconstruction in electric circuit networks, Phys. Rev. B 99, 161114(R) (2019).
  * (18) S. M. Rafi-Ul-Islam, Z. B. Siu, and M. B. A. Jalil, Topoelectrical circuit realization of a Weyl semimetal heterojunction, Commun. Phys. 3, 72 (2020).
  * (19) J. Dong, V. Juričić, and B. Roy, Topolectric circuits: Theory and construction, Phys. Rev. Research 3, 023056 (2021).
  * (20) Z.-Q. Zhang, B.-L. Wu, J. Song, and H. Jiang, Topological Anderson insulator in electric circuits, Phys. Rev. B 100, 184202 (2019).
  * (21) E. Zhao, Topological circuits of inductors and capacitors, Ann. Phys. 399, 289 (2018).
  * (22) Y. Nakata, T. Okada, T. Nakanishi, and M. Kitano, Circuit model for hybridization modes in metamaterials and its analogy to the quantum tight binding model, Phys. Status Solidi B 249, 2293 (2012).
  * (23) M. Ezawa, Higher-order topological electric circuits and topological corner resonance on the breathing kagome and pyrochlore lattices, Phys. Rev. B 98, 201402(R) (2018).
  * (24) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, and T. Neupert, Topolectrical-circuit realization of topological corner modes, Nat. Phys. 14, 925 (2018).
  * (25) R. Haenel, T. Branch, and M. Franz, Chern insulators for electromagnetic waves in RLC networks, Phys. Rev. B 99, 235110 (2019).
  * (26) W. Zhu, Y. Long, H. Chen, and J. Ren, Quantum valley Hall effects and spin-valley locking in topological Kane-Mele circuit networks, Phys. Rev. B 99, 115410 (2019).
  * (27) L. Song, H. Yang, Y. Cao, and P. Yan, Realization of the Square-Root Higher-Order Topological Insulator in Electric Circuits, Nano Lett. 20(10), 7566 (2020).
  * (28) S. M. Rafi-Ul-Islam , Z. B. Siu, and M. B. A. Jalil, Topological phases with higher winding numbers in nonreciprocal one-dimensional topolectrical circuits, Phys. Rev. B 103, 035420 (2021).
  * (29) F. Y. Wu, Theory of resistor networks: the two-point resistance, J. Phys. A: Math. Gen. 37, 6653 (2004).