op-symbol = d, *-op-left = true, op-order-sep = 0 mu,

# Analysis of autocorrelation times in Neural Markov Chain Monte Carlo
simulations

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

We provide a deepened study of autocorrelations in Neural Markov Chain Monte
Carlo simulations, a version of the traditional Metropolis algorithm which
employs neural networks to provide independent proposals. We illustrate our
ideas using the two-dimensional Ising model. We propose several estimates of
autocorrelation times, some inspired by analytical results derived for the
Metropolized Independent Sampler, which we compare and study as a function of
inverse temperature $\beta$. Based on that we propose an alternative loss
function and study its impact on the autocorelation times. Furthermore, we
investigate the impact of imposing system symmetries ($Z_{2}$ and/or
translational) in the neural network training process on the autocorrelation
times. Eventually, we propose a scheme which incorporates partial heat-bath
updates. The impact of the above enhancements is discussed for a $16\times 16$
spin system. The summary of our findings may serve as a guidance to the
implementation of Neural Markov Chain Monte Carlo simulations of more
complicated models.

##  1 Introduction

The original idea behind Monte Carlo simulations proposed by Stanislaw Ulam
and Enrico Fermi and later expressed in a form of an algorithm applied to
study a simple classical statistical mechanics problem by Metropolis et al.
[1] turned out to be a very powerful approach to tackle a large variety of
problems in all computational sciences (see for example [2]). In case of
complicated probability distributions most formulations resort to the
construction of an associated Markov chain of consecutive proposals. The
statistical uncertainty of any outcome of Monte Carlo simulation depends
directly on the number of statistically independent configurations used to
estimate it. Hence, the effectiveness of simulation algorithms is measured by
autocorrelation time which quantify how many configurations are produced by
the algorithm before a new, statistically independent configuration appears.
Increasing autocorrelation times, a phenomenon called critical slowing down,
is usually the main factor which limits the statistical precision of outputs.
In the context of field theory simulations in elementary particle physics
several proposals were advanced in order to alleviate that problem, in
particular metadynamics [3, 4], instanton updates [5] or multiscale
thermalization [6].

The recent great interest in machine learning techniques has also provided new
ideas in the domain of Monte Carlo simulations which aim at reducing the
autocorrelation times. The ability of artificial neural networks to
approximate a very wide class of probability distributions was used in Ref.[7]
to propose a variational estimate of free energy of statistical systems.
Subsequently the idea was extended and used as a mechanism of providing
uncorrelated proposals in a Monte Carlo simulation in Ref.[8]. The Authors
call the resulting algorithm Neural Markov Chain Monte Carlo (NMCMC).

Since the proposal of the NMCMC algorithm a few studies discussed its
effectiveness. In particular, the Authors of Ref. [9] applied a similar method
to spin glass models, whereas in Ref. [10] an alternative update mechanism was
proposed. In the present study, we concentrate on the classical test bench
provided by the two-dimensional Ising model in order to stay in a simple and
well understood physical situation and focus our attention on the algorithmic
effectiveness of the approach. We quantify the latter using several possible
measures of autocorrelation times and compare them with the original
Metropolis algorithm [1] and with the cluster algorithm proposed by Wolff [11,
12]. The Ising model offers also exact solution for the partition function of
the finite-size system [13] from which the free energy can be obtained and
which provides measures of how well the neural network approximates the
desired probability distribution.

In the following we start in Section 1.1 by briefly defining the model, in
Section 1.2 we present the Neural Markov Chain Monte Carlo (NMCMC) algorithm
using Autoregressive Networks (AN)[8]. In Section 2 we summarize the known
analytical results for the autocorrelation time in a simulation using the
Metropolized Independent Sampler (MIS) derived by Liu in Ref. [14] and discuss
them in the present setup. Subsequently, in Section 3 we define several
possible estimators for the autocorrelation time. In Section 4 we collect two
loss functions used in the present context in the Literature and we supplement
them with a new loss function inspired by the analytic insight discussed in
Section 2. The results of our numerical experiments are divided into two
parts: Sections 5 and 6. In the former part we use a small system size which
enables us to estimate all analytical quantities of interest and provide their
comparisons with the approximated counterparts defined on a small number of
configurations contained in a batch. In the latter part we use a larger system
size in order to check the algorithmic effectiveness for more realistic system
sizes. Finally, in the remaining Section 7 we describe two techniques which
improve the quality of the training of the neural network. First, we exploit
the symmetries. In Section 7.1 we monitor the effect of enforcing the $Z_{2}$
and translational symmetries during the training process. Second, in Section
7.2, we use the locality of nearest-neighbour interactions of the Hamiltonian
of the Ising model to reduce the number of relevant spin that have to be
modelled by the neural network. By using so-called chessboard partitioning
half of the spins has to be generated by the neural network while the
remaining half can be obtained exactly from the local Boltzmann probability
distribution providing the heath-bath type of update. We provide a summary of
our results and an outlook in Section 8.

###  1.1 Ising model

We consider the two-dimensional ferromagnetic Ising model with $N=L^{2}$ spins
which are located on the $L\\!\\!\times\\!\\!L$ square lattice. We will always
compute observables on configurations coming from a Markov chain, hence
labelled by a ”Monte Carlo” time. We denote the $k$-th configuration in the
chain by $\mathbf{s}_{k}$. The spins within a single configuration are
labelled by a linear upper index $i$; thus the $i$-th spin of $k$-th
configuration is $s_{k}^{i}$. Each state of the system is thus a configuration
of spins $\mathbf{s}_{k}=(s_{k}^{1},\ldots s_{k}^{N})$, and $s_{k}^{i}=\pm 1$.
For a configuration of spins $\mathbf{s}_{k}$ the energy of that configuration
is given by the Hamiltonian

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

where the sum runs over all neighbouring pairs of lattice sites. In the
following, we set the coupling constant $J=1$ without loosing generality. We
assume periodic boundary conditions for the lattice.

The corresponding Boltzmann distribution is then given by

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

with partition function $Z=\sum_{\mathbf{s}}\exp(-\beta H(\mathbf{s}))$, where
the sum is performed over all $M=2^{N}$ states of the system.

It is a well known fact that, in the infinite volume limit, the model
undergoes second-order phase transition at
$\beta_{c}=\frac{1}{2}\ln(1+\sqrt{2})\approx 0.441$. The corresponding order
parameter is the absolute magnetization per spin

$m=\left\langle{\frac{1}{N}\left|\sum_{i=1}^{N}\mathbf{s}^{i}\right|}\right\rangle_{p}=\sum_{k=1}^{M}p(\mathbf{s}_{k})\left(\frac{1}{N}\left|\sum_{i=1}^{N}\mathbf{s}^{i}_{k}\right|\right)\approx\frac{1}{N_{batch}}\sum_{k=1}^{N_{batch}}\frac{1}{N}\left|\sum_{i=1}^{N}\mathbf{s}^{i}_{k}\right|,$
(3)

where the last expression is an average over $N_{batch}\gg 1$ configurations
sampled from distribution $p$. In a finite volume the system exhibits a
crossover at $\beta_{c}(L)$, i.e. the value of $\beta$ at given $L$ where
$|\frac{\partial m}{\partial\beta}|$ is maximal. We assume that the finite
volume system is in an ordered state when $\beta>\beta_{c}(L)$ and in a
disordered state for $\beta<\beta_{c}(L)$.

###  1.2 Neural Markov Chain Monte Carlo (NMCMC)

In the approach discussed in this paper, we use an artificial neural network
to model the probability $p(\mathbf{s})$. In order to do so we must first
recast $p(\mathbf{s})$ in the form that will ensure proper normalisation and
allow for easy sampling from the learned distribution. To this end we rewrite
$p(\mathbf{s}_{k})$ as a product of conditional probabilities of consecutive
spins, i.e.

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

Note that factorization of the total probability of the configuration into $N$
conditional probabilities requires a unique ordering of all spins on the
lattice. However, the precise form the ordering does not matter if it is
maintained during the entire simulation. For definiteness, in our
implementation we use a lexicographic ordering on the two-dimensional lattice,
the $x$ coordinate being faster than the $y$ coordinate.

Now we can use an artificial neural network to model the conditional
probabilities $p(s^{i}|\cdot)$. In the simplest realization, the neural
network will have $N$ input neurons and $N$ output neurons. The value of the
output of the $i$-th neuron is interpreted as

$q_{\theta}(s^{i}=1|s^{1},s^{2},\dots,s^{i-1})$ (5)

In order to ensure that the conditional probabilities depend only on the
values of the preceding spins a particular architecture of artificial neural
networks is used, the so-called autoregressive networks. They are build out of
a class of dense or convolutional layers where the neural connections are only
allowed to the neurons associated to the preceding spins. The details of such
a construction can be found in Ref.[7].

It is the purpose of the training steps to tune the weights of the neural
network, collectively denoted as $\theta$ in such a way as to approximate
$p(\mathbf{s})$ as close as possible. This is accomplished by looking for the
minimum of the loss function, usually the Kullback–Leibler divergence which
quantifies the difference of two probability distributions. In the following
we also discuss other possible choices and their interpretation.

We can easily generate configurations from the $q_{\theta}(\mathbf{s})$
distribution by generating subsequent spins from conditional probabilities
$q_{\theta}(s^{i}|s^{1},\ldots,s^{i-1})$. This requires $N$ sequential calls
to $q_{\theta}(s|\cdot)$ functions i.e. making $N$ forward passes through the
neural network. Please note that this procedure does not rely in any place on
any preceding spin configuration, hence it provides a mechanism for generating
uncorrelated configurations from the $q_{\theta}$ distribution avoiding the
construction of a Markov chain.

In the end we are interested in spin configurations generated with the target,
Boltzmann probability distribution. Having a neural network generating states
according to the distribution $q_{\theta}$ we can correct the difference
between $p(\mathbf{s})$ and $q_{\theta}(\mathbf{s})$ within an accept/reject
step with the accept probability for the state $\mathbf{s}_{k+1}$,

$\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).$
(6)

This is know as metropolized independent sampler [14, 15, 16, 17].

The fraction of the accepted propositions is another measure of the difference
between $q_{\theta}$ and $p$. When $q_{\theta}=p$ acceptance ratio is 1. When
two distributions are not identical some propositions can be rejected what
leads to repetitions in the Markov chain. Therefore, even though propositions
are independent from each other, a nonzero autocorrelation is possible. The
range of those correlations is determined by the eigenvalues and eigenvectors
of the transition matrix of the Markov-Chain which are the subject of the next
section.

##  2 Analytic results for the eigenvalues of transition matrix

For the sake of completeness we include some results from Ref.[14] and we use
this opportunity to set up the notation. We consider two probability
distributions: $q(\mathbf{s})$ and $p(\mathbf{s})$, defined on some set $X$.
Distribution $p(\mathbf{s})$ is the target probability distributions. Using
the metropolised independent sampler as described in the previous section to
sample from $p$ we generate a configuration $\mathbf{s}_{k+1}$ from $q$ and
accept it with probability

$\mathcal{P}(\mathbf{s}_{k}\rightarrow\mathbf{s}_{k+1})=\min\left\\{1,\frac{w(\mathbf{s}_{k+1})}{w(\mathbf{s}_{k})}\right\\},$
(7)

where the importance ratio $w(\mathbf{s})$ is defined for each configuration
$\mathbf{s}\in X$ by

$w(\mathbf{s})=\frac{p(\mathbf{s})}{q(\mathbf{s})}.$ (8)

Since $p(\mathbf{s}),q(\mathbf{s})>0$ it follows that
$0<w(\mathbf{s})<\infty$111In principle $q(s)$ could become zero for some $s$,
but then it would mean that this configuration would never be proposed as the
trial configuration thus avoiding the problem of infinite $w(\mathbf{s})$.

The main result of Liu presented in Ref.[14] provides analytic formulae for
the eigenvalues of the transition matrix of the resulting Markov chain in
terms of the importance ratios. For a system with a finite number $M$ of
elements in $X$ one can order the importance ratios $w(\mathbf{s})$ in an
decreasing order (countable and continuous sets $X$ were also discussed in
[14]), so we adopt the labeling (which changes if $p$ or $q$ is modified) of
configurations $\mathbf{s}_{i}$ such that the corresponding importance ratios
satisfy

$w(\mathbf{s}_{1})\geq w(\mathbf{s}_{2})\geq\dots\geq w(\mathbf{s}_{M}).$ (9)

Direct diagonalization of the transition matrix of the Markov chain gives its
eigenvalues $\lambda_{k}$, where

$\displaystyle\lambda_{k}=\left\\{\begin{array}[]{cc}1&\textrm{for }k=0,\\\
\sum\limits_{i=k}^{M}q(\mathbf{s}_{i})\left(1-\frac{w(\mathbf{s}_{i})}{w(\mathbf{s}_{k})}\right)&\textrm{for
}0<k\leq M-1\end{array}\right.$ (12)

The eigenvalue $\lambda_{0}=1$ corresponds to the stationary state and all
remaining eigenvalues $\lambda_{k>0}$ are nonnegative and smaller than 1. Due
to the introduced ordering of importance ratios (9) and the fact that all $w$
are positive, it follows that
$1>\lambda_{1}>\lambda_{2}>\ldots>\lambda_{M-1}>0$.

In the case of the Neural Markov Chain Monte Carlo algorithm, the sampling
distribution $q(\mathbf{s})$ is encoded in the neural network and therefore
depends on all the weights $\theta$ of the network,
$q(\mathbf{s})=q_{\theta}(\mathbf{s})$. The aim of the training process is to
ensure that $q_{\theta}(\mathbf{s})$ is as close to $p(\mathbf{s})$ as
possible thus driving the importance ratios closer to 1 and hence decreasing
the values of $\lambda_{k>0}$. This in turn translates into the reduction of
the autocorrelation times as explained in more details in the next section.

|
---|---
Figure 1: Histograms of $M-1$ eigenvalues of transition matrix,
$\lambda_{k>0}$, for system $4\times 4$. Training was performed using the
$D_{\textrm{KL}}$ loss function. Left figure is for initial state of network,
right is for fully trained network. Green line denotes $M-1=2^{16}-1=65535$
value.

In Fig. 1 we show histograms of all eigenvalues of transition matrix smaller
than 1 – given by eq. (12) for Ising model of size $4\times 4$ (so
$2^{16}-1=65535$ eigenvalues in total). Temperature was chosen as $\beta=0.6$
and the $D_{\textrm{KL}}$ loss function was employed (see section 4). Interval
$(0,1)$ was divided into 100 bins.

The left histogram shows the situation when the network is at the initial
state (1 epoch), right panel shows the histogram when the neural network is
trained. At the beginning of the simulation (1 epoch) the eigenvalues are
distributed randomly with fluctuations coming from the random initial network
weights. During the training process the eigenvalues decrease, as expected,
since $q_{\theta}$ approaches $p$ distribution. The largest eigenvalue,
$\lambda_{1}$ is around 0.75 when the network is fully trained.

##  3 Autocorrelation times

In order to quantitatively discuss the algorithmic effectiveness we will
estimate and compare the associated autocorrelation times. In practice there
exist several ways to define the relevant autocorrelation time. On one hand,
one can focus on the transition matrix of the Markov chain and consider the
analytically known or approximated autocorrelation time extracted from the
eigenvalues of that matrix (see Section 2). On the other hand, one can take
the perspective of a particular observable of interest and estimate the
autocorrelation time of that particular observable. In certain limiting
situations all mentioned approaches should provide comparable estimates of the
same, say longest, autocorrelation time. However, in the case of finite Markov
chains and finite statistical precision, the autocorrelation times extracted
along the three ways may vary by several orders of magnitude. We explain this
effect in Section 5, meantime we provide all needed definitions below. The
autocorrelation function is defined as

$\Gamma_{\mathcal{O}}(t)=\frac{1}{\operatorname{var}\,\mathcal{O}}E\left[\left(\mathcal{O}(\mathbf{s}_{i})-\bar{\mathcal{O}}\right)\left(\mathcal{O}(\mathbf{s}_{i+t})-\bar{\mathcal{O}}\right)\right],$
(13)

where the expectation value is calculated with respect to all the realisations
of the Markov chain. The $\bar{\mathcal{O}}$ is the average of the observable
$\mathcal{O}$ with respect to the equilibrium distribution of the Markov
chain. Normalisation is such that $\Gamma(0)=1$. This function can be solely
determined from the properties of the eigenvectors and eigenvalues of the
Markov chain transition matrix. In general autocorrelation function
$\Gamma(t)$ for any observable $\mathcal{O}$ will be of the form:

$\Gamma_{\mathcal{O}}(t)=\sum_{k>0}c_{\mathcal{O},k}^{2}\,\lambda_{k}^{t}=\sum_{k>0}c_{\mathcal{O},k}^{2}\,e^{-\frac{t}{\tau_{k}}},\quad\text{with
}\tau_{k}=-(\log\lambda_{k})^{-1}.$ (14)

The coefficients $c_{\mathcal{O},k}$ can be calculated using the eigenvectors
of transition matrix. One expects that for large separations $t$ the
autocorrelation function decays exponentially:

$\Gamma_{\mathcal{O}}(t)\sim e^{-t/{\tau_{L}^{\mathcal{O}}}}\,\textrm{ for
}t\rightarrow\infty.$ (15)

The $\tau_{L}^{\mathcal{O}}$ will be given by the largest eigenvalue
$\lambda_{k}$ for which $c_{\mathcal{O},k}>0$, and we can define

$\tau_{L}=\sup_{O}\tau_{L}^{\mathcal{O}}=-(\log\lambda_{1})^{-1},$ (16)

where $\sup_{O}$ is taken over all possible observables.

In the extreme case when $w(\mathbf{s})=1$ for all $\mathbf{s}$, i.e. $q=p$,
it follows that $\lambda_{k>0}=0$ and therefore $\tau_{L}=0$. Note that we do
not refer to any particular observable in defining $\tau_{L}$. It is a
property of the Markov chain itself. Different observables will couple to the
different modes of the Markov chain and will exhibit different leading
autocorrelation times, which however will be always smaller or equal to
$\tau_{L}$.

The most important measure of autocorrelations is the integrated
autocorrelation time which enters as a factor into the definition of the
statistical uncertainty of any value estimated from the Markov chain. For the
observable $\mathcal{O}$ it is defined as [18]:

$\tau^{O,int}=1+2\sum_{i=1}^{\infty}\Gamma_{\mathcal{O}}(t).$ (17)

In order to operate at the level of integrated autocorrelation times we can
define ”integrated version” of the $\tau_{L}^{\mathcal{O}}$:

$\tau_{L}^{\mathcal{O},int}\equiv
1+2\sum_{t=1}^{\infty}e^{-\frac{t}{\tau_{L}^{\mathcal{O}}}}=\frac{1+e^{-{1}/{\tau_{L}^{\mathcal{O}}}}}{1-e^{-{1}/{\tau_{L}^{\mathcal{O}}}}}$
(18)

and similarly for $\tau_{L}^{int}$

$\tau_{L}^{int}=\frac{1+e^{-{1}/{\tau_{L}}}}{1-e^{-{1}/{\tau_{L}}}}.$ (19)

Note that among those three theoretical autocorrelation times only
$\tau_{L}^{int}$ can be calculated without knowing $c_{\mathcal{O},k}$
coefficients. What is more, this is possible only for very small systems where
all states can be constructed explicitly. Calculating $c_{\mathcal{O},k}$ is
beyond scope of this paper. Instead, we define below other estimators of
integrated autocorrelation times which are calculable for systems of any size.

###  3.1 Estimators of autocorrelation times

In practice we estimate $\Gamma_{\mathcal{O}}(t)$ from the finite sample:

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

where $\bar{\mathcal{O}}$ and $\operatorname{var}\,\mathcal{O}$ are now mean
and variance calculated using $N_{sample}$ states. In the following we mark
all the quantities that are estimated from finite sample with a hat sign. The
autocorrelation function $\Gamma_{\mathcal{O}}(t)$ is always nonnegative but
its estimator $\hat{\Gamma}_{\mathcal{O}}(t)$ can become negative for some
$t_{max}$ due to finite fluctuations. For the purpose of these studies we
consider autocorrelation function for separations smaller than $t_{max}$:
$\hat{\Gamma}_{\mathcal{O}}(t),\ t<t_{max}$. This immediately gives us the
estimator of integrated time Eq. (17):

$\hat{\tau}^{\mathcal{O},int}=1+2\sum_{i=1}^{t_{max}}\hat{\Gamma}_{\mathcal{O}}(t).$
(21)

To estimate $\tau_{L}^{\mathcal{O},int}$ we fit the exponent to the
autocorrelation function $\hat{\Gamma}$ at large $t$. In order to take into
account also subleading autocorrelation time which can contribute at moderate
$t$, we use two exponent function of the form:

$a_{1}e^{-t/\tau_{1}}+a_{2}e^{-t/\tau_{2}}.$ (22)

We have checked that introducing more terms in Eq. (22) does not change the
value of largest $\tau_{i}$. We thus define
$\hat{\tau}^{\mathcal{O}}_{F}=\max(\tau_{1},\tau_{2})$ which describes the
behaviour of $\hat{\Gamma}_{\mathcal{O}}(t)$ at large $t$ (but $t<t_{max}$).
The integrated version of $\hat{\tau}^{\mathcal{O}}_{F}$ is:

$\hat{\tau}^{\mathcal{O},int}_{F}=\frac{1+e^{-{1}/{\hat{\tau}^{\mathcal{O}}_{F}}}}{1-e^{-{1}/{\hat{\tau}^{\mathcal{O}}_{F}}}}.$
(23)

The quantity $\hat{\tau}^{\mathcal{O},int}_{F}$ given by this definition is an
estimator of $\tau_{L}^{\mathcal{O},int}$ defined in Eq. (18).

Also definition Eq. (16) of $\tau_{L}$ has an estimator that is readily
available and can be monitored during the training process of the neural
network. The exact value of $\lambda_{1}$ is not available for most systems,
but can be approximated by $\hat{\lambda}_{1}$ calculated from a finite sample
of configurations (also called batch). Hence, we define the ”batch version” of
largest eigenvalue $\lambda_{1}$,

$\hat{\lambda}_{1}\equiv\frac{1}{N_{\textrm{batch}}}\sum_{k=1}^{N_{\textrm{batch}}}\Big{(}1-\frac{w({\mathbf{s}_{k}})}{w_{\textrm{max
on batch}}}\Big{)}\Big{|}_{\textrm{with }\mathbf{s}_{k}\textrm{ sampled from
}q_{\theta}},$ (24)

and correspondingly the ”batch version” of leading autocorrelation time
$\tau_{L}$:

$\hat{\tau}_{B}=-(\log\hat{\lambda}_{1})^{-1},$ (25)

which is an analogue of Eq. (16). Integrated version of this $\hat{\tau}_{B}$
is:

$\hat{\tau}_{B}^{int}\equiv\frac{1+e^{-{1}/{\hat{\tau}_{B}}}}{1-e^{-{1}/{\hat{\tau}_{B}}}}.$
(26)

To conclude, we have defined here $\hat{\tau}^{\mathcal{O},int}$,
$\hat{\tau}^{\mathcal{O},int}_{F}$, $\hat{\tau}_{B}^{int}$ which are
estimators of three versions of theoretical integrated autocorrelation times:
$\tau^{O,int}$, $\tau_{L}^{\mathcal{O},int}$ and $\tau_{L}^{int}$. In the
following we choose the energy $H$, Eq. (1), as observable $\mathcal{O}$. We
have checked that the absolute value of the magnetization gives very similar
results. We estimate $\hat{\Gamma}_{\mathcal{O}}(t)$ at the fixed stage of the
training process and we generate a long, $\mathcal{O}(10^{5})$, sequence of
states. Subsequently we apply the condition Eq. (6) for consecutive states in
this sequence and obtain a Markov chain whose stationary state is the
Boltzmann distribution.

###  3.2 Some relations among the different definitions

Let us now comment on the hierarchy of (integrated) autocorrelation times we
have just defined, $\tau^{O,int}$, $\tau_{L}^{\mathcal{O},int}$ and
$\tau_{L}^{int}$. By construction the largest of them will be
$\tau_{L}^{int}$. As mentioned earlier, it is the property of the Markov chain
not associated to any particular observable. Hence,
$\tau_{L}^{\mathcal{O},int}$ and $\tau^{O,int}$ are smaller than or equal to
$\tau_{L}^{int}$. In principle $\tau^{O,int}$ is equal to $\tau_{L}^{int}$
only if the observable $\mathcal{O}$ couples exclusively to the slowest mode
of the Markov chain. One can show that once the observable has a non-zero
coupling to another mode, $\tau^{O,int}$ will be strictly smaller than
$\tau_{L}^{int}$. Also by definition, $\tau^{O,int}$, being the integral of
the autocorrelation function, is smaller than $\tau_{L}^{\mathcal{O},int}$,
except for the case where the observable couples exactly to one of the modes
of the Markov chain. In that situation we expect that
$\tau^{O,int}=\tau_{L}^{\mathcal{O},int}$.

As for the $\hat{\tau}^{\mathcal{O},int}$, $\hat{\tau}^{\mathcal{O},int}_{F}$,
$\hat{\tau}_{B}^{int}$ estimators their hierarchy cannot be theoretically
established as they are all calculated in a very different way. However we
could expect that the ordering presented above is, at least approximately,
preserved leading to
$\hat{\tau}_{B}^{int}\geq\hat{\tau}^{\mathcal{O},int}_{F}\geq\hat{\tau}^{\mathcal{O},int}$.
We present the measurements of the above integrated autocorrelation times on a
$4\times 4$ and $16\times 16$ lattice for the Ising model in the following
sections and find that while
$\hat{\tau}^{\mathcal{O},int}_{F}\geq\hat{\tau}^{\mathcal{O},int}$,
$\hat{\tau}_{B}^{int}$ can be smaller than $\hat{\tau}^{\mathcal{O},int}_{F}$.
We would also expect that all autocorrelation times are smaller than
$\tau_{L}^{int}$ and this is indeed the case.

##  4 Neural network loss functions

We explained in the previous section that the difference between the target
probability distribution $p(\mathbf{s})$ and the sampling probability
distribution $q_{\theta}(\mathbf{s})$ directly affects the autocorrelation
times via Eq.(12). Several measures can be used to quantify the distance
between two distributions, in the context of machine learning the most popular
being the Kullback-Leibler divergence. Of course, in the ideal case, finding
the global zero of any of the correctly defined measures will lead to finding
such $\theta$ weights that $q_{\theta}(\mathbf{s})=p(\mathbf{s})$. However,
for a finite size of the neural network and thus at a limited ability of the
neural networks to model any probability function as well as for a limited
number of training epochs, different measures may potentially be more
efficient than others as they may lead to $q_{\theta}(\mathbf{s})$ that gives
smaller autocorrelation times. Below we provide three definitions of loss
functions which we will apply to the case of two dimensional Ising model
following sections.

###  4.1 Kullback–Leibler (KL) divergence

The most commonly used in the literature loss function is defined as:

$\displaystyle D_{\textrm{KL}}(q_{\theta}|p)$
$\displaystyle=\sum_{k=1}^{M}q_{\theta}(\mathbf{s}_{k})\,\log\left(\frac{q_{\theta}(\mathbf{s}_{k})}{p(\mathbf{s}_{k})}\right).$
(27)

Substituting Eq.(2) into the last equation one obtains:

$\displaystyle D_{\textrm{KL}}(q_{\theta}|p)$ $\displaystyle=\beta(F_{q}-F),$
(28)

where $F=-\frac{1}{\beta}\log Z$ is a true free energy and

$F_{q}=\frac{1}{\beta}\sum_{k=1}^{M}q_{\theta}(\mathbf{s}_{k})\left[\beta
H(\mathbf{s}_{k})+\log q_{\theta}(\mathbf{s}_{k})\right]$ (29)

is variational free energy. Since $D_{\textrm{KL}}$ is a non-negative quantity
the minimalization of KL divergence is equivalent to minimization of the
variational free energy.

As far as the importance ratios are concerned we can easily identify their
presence in Eq.(27), since we can write

$D_{\textrm{KL}}(q_{\theta}|p)=-\sum_{k=1}^{M}q_{\theta}(\mathbf{s}_{k})\,\log
w(\mathbf{s}_{k})=-\left\langle{\log w}\right\rangle_{q_{\theta}}.$ (30)

It should be clear from this reformulation that during the training process
the neural network will try to improve $q_{\theta}$ for configurations for
which the importance ratio is the smallest, i.e. for which
$q_{\theta}(\mathbf{s}_{k})$ is drastically overestimated with respect to
$p(\mathbf{s}_{k})$.

Calculating the gradient of the loss function

$\diff{D_{\textrm{KL}}(q_{\theta}|p)}{\theta}=-\diff*{\left\langle{\log
w}\right\rangle_{q_{\theta}}}{\theta}$ (31)

requires certain caution. The derivative in Eq. (31) acts not only on the
expression being averaged but also on the distribution used for averaging.
More generally

$\diff{}{\theta}\Big{\langle}S(\mathbf{s}|\theta)\Big{\rangle}_{q_{\theta}}=\diff*{\sum_{\mathbf{s}}{q_{\theta}}(\mathbf{s})S(\mathbf{s}|\theta)}{\theta}=\left\langle{\diff{\log{q_{\theta}}}{\theta}S}\right\rangle_{q_{\theta}}+\left\langle{\diff{S}{\theta}}\right\rangle_{q_{\theta}}.$
(32)

In the case of KL divergence, $S=\log q_{\theta}-\log p$, and

$\left\langle{\diff{S}{\theta}}\right\rangle_{q_{\theta}}=\diff*{\sum_{s}{q_{\theta}}(s)}{\theta}=0.$
(33)

So finally

$\diff{D_{\textrm{KL}}(q_{\theta}|p)}{\theta}=\left\langle{\diff{\log
q_{\theta}}{\theta}(\log q_{\theta}-\log p)}\right\rangle_{q_{\theta}}.$ (34)

As was noticed in Ref. [7], Eq. (34) can be interpreted as a REINFORCE
algorithm [19] with score function $-S$.

###  4.2 Largest Eigenvalue of transition matrix $\lambda_{1}$

An alternative way to measure the difference between $q_{\theta}$ and $p$ is
to use directly the expression for the largest eigenvalue of the associated
Markov chain transition matrix provided by Eq. (12)

$D_{\lambda_{1}}(q_{\theta}|p)\equiv\lambda_{1}=\sum_{k=1}^{M}q_{\theta}(\mathbf{s}_{k})\left(1-\frac{w(\mathbf{s}_{k})}{w(\mathbf{s}_{1})}\right)=\left\langle{1-\frac{w}{w(\mathbf{s}_{1})}}\right\rangle_{q_{\theta}}.$
(35)

Because $\left\langle{w}\right\rangle_{q_{\theta}}=1$

$D_{\lambda_{1}}(q_{\theta}|p)=1-\frac{1}{w(\mathbf{s}_{1})}.$ (36)

We cannot however use this form without knowing the the overall normalization
factor, partition function $Z$. In (35) $p$ enters only as a ratio and so the
partition function $Z$ is not needed to calculate the eigenvalue. As already
stated, $\lambda_{1}\geq 0$ and it vanishes only when $q_{\theta}\equiv p$.

Although this loss function has the same minimum as the KL divergence, the
contributions coming from particular configurations are different, and hence
during minimization $q_{\theta}$ will be brought closer to $p$ for a different
subset of configurations first.

Using Eq. (32) one can obtain the gradient of $D_{\lambda_{1}}$,

$\displaystyle\diff{D_{\lambda_{1}}(q_{\theta}|p)}{\theta}=$
$\displaystyle\left\langle{\diff{\log
q_{\theta}}{\theta}}\right\rangle_{q_{\theta}}-\left\langle{\frac{w}{w(\mathbf{s}_{1})}\diff{\log
q_{\theta}(\mathbf{s}_{1})}{\theta}}\right\rangle_{q_{\theta}}.$ (37)

Alternatively, one can follow the REINFORCE algorithm and use the following
quantity in backward propagation,

$\left\langle{\diff{\log
q_{\theta}}{\theta}\left(1-\frac{w}{w(\mathbf{s}_{1})}\right)}\right\rangle_{q_{\theta}},$
(38)

instead of the derivative Eq. (37). Note that, in contrast to
$D_{\textrm{KL}}$, these two approaches are different. Numerical studies
showed that the approach given by Eq. (38) results in a much more efficient
training for the neural network. The comparison of these two approaches and
their efficiencies is beyond scope of this paper. In what follows we show
results only for the REINFORCE approach.

###  4.3 Effective Sample Size (ESS)

Following Ref.[14, 20] we define yet another possible loss function, which for
the sake of completeness we compare with the two previously defined,

$ESS=\frac{\left\langle{w}\right\rangle_{q_{\theta}}^{2}}{\left\langle{w^{2}}\right\rangle_{q_{\theta}}}.$
(39)

The quantity is normalized such that $ESS\in[0,1]$ and for perfectly matched
distributiond one has $ESS=1$. So we define the loss function as $1-ESS$

$D_{\textrm{ESS}}(q_{\theta}|p)=1-\frac{\left\langle{w}\right\rangle_{q_{\theta}}^{2}}{\left\langle{w^{2}}\right\rangle_{q_{\theta}}}=\frac{\left\langle{w^{2}}\right\rangle_{q_{\theta}}-\left\langle{w}\right\rangle_{q_{\theta}}^{2}}{\left\langle{w^{2}}\right\rangle_{q_{\theta}}}=\frac{\operatorname{var}_{q_{\theta}}{w}}{\left\langle{w^{2}}\right\rangle_{q_{\theta}}}.$
(40)

Again, the partition function $Z$ is not needed to calculate $ESS$.

This loss function cannot be expressed as a weigthed average so the REINFORCE
algorithm is not available. We differentiate Eq. (40) obtaining:

$\diff{D_{\textrm{ESS}}(q_{\theta}|p)}{\theta}=-\frac{\left\langle{w}\right\rangle_{q_{\theta}}^{2}}{\left\langle{w^{2}}\right\rangle_{q_{\theta}}^{2}}\left\langle{\diff{\log{q_{\theta}}}{\theta}\;w^{2}}\right\rangle_{q_{\theta}}.$
(41)

###  4.4 Batch versions of the loss functions

For all but very small systems, the sum over all configurations
$\sum_{i=k}^{M}(\ldots)$ is prohibitively large to perform exactly. In the
training process one disposes of a set of configurations $\mathbf{s}$ drawn
from the current probability distribution $q_{\theta}(\mathbf{s})$. Hence, all
the three loss functions defined above can be expressed as batch averages. To
this goal we do the following replacements in Eqs. (30), (35) and (40):

$\left\langle{\ldots}\right\rangle_{q_{\theta}}=\sum_{i=1}^{M}q_{\theta}(\mathbf{s}_{i})(\ldots)\rightarrow\frac{1}{N_{batch}}\sum_{i=1}^{N_{batch}}(\ldots).$
(42)

For the KL divergence we thus have,

$D_{\textrm{KL}}(q_{\theta}|p)\approx\hat{D}_{\textrm{KL}}(q_{\theta}|p)\equiv-\frac{1}{N_{\textrm{batch}}}\sum_{k=1}^{N_{\textrm{batch}}}\log
w(\mathbf{s}_{k})\Big{|}_{\textrm{with }\mathbf{s}_{k}\textrm{ sampled from
}q_{\theta}}=-\langle\log w\rangle_{\textrm{batch average}},$ (43)

where we used the hat symbol to denote quantity calculated from the batch of
configurations. In the case of $D_{\lambda_{1}}(q_{\theta}|p)$ we also need to
modify the value of $w_{1}$ since in practice its value may be unknown. We
thus say

$\begin{split}D_{\lambda_{1}}(q_{\theta}|p)\approx\hat{D}_{\lambda_{1}}(q_{\theta}|p)&\equiv\frac{1}{N_{\textrm{batch}}}\sum_{k=1}^{N_{\textrm{batch}}}\Big{(}1-\frac{w({\mathbf{s}_{k}})}{w_{\textrm{max
on batch}}}\Big{)}\Big{|}_{\textrm{with }\mathbf{s}_{k}\textrm{ sampled from
}q_{\theta}}\\\ &=\langle 1-\frac{w}{w_{\textrm{max on
batch}}}\rangle_{\textrm{batch average}}\end{split}$ (44)

It should be clear that $w(\mathbf{s}_{1})\geq w_{\textrm{max on batch}}$.
This however does not necessarily entails that
$\hat{D}_{\lambda_{1}}(q_{\theta}|p)$ is smaller than
$D_{\lambda_{1}}(q_{\theta}|p)$ as the relation Eq.(36) need not hold for the
finite sample.

##  5 Small lattice

In Section 2 we recalled analytical results concerning NMCMC algorithms, in
particular explicit expressions for the eigenvalues of the transition matrix.
The evaluation of these formulae is only possible when the system size is
small and one can visit all configurations $\mathbf{s}\in X$. The aim of this
section is to study a system small enough, we choose the lattice of size
$4\times 4$, and investigate the differences between the exact analytic
results and their approximated batch versions. For systems of physically
reasonable sizes only the latter are available and therefore it is important
to understand their properties.

###  5.1 Neural network architecture

In this work we employ the simplest implementation of an autoregressive neural
network with dense layers. We use one hidden layer and the width of all layers
is equal to the number of spins on the lattice. As activation functions we use
PReLU in between layers and sigmoid for output layer. Adam optimizer is used.
The $\beta$-annealing is used for training, except for Figs. 1, 2, 3. We
demonstrate in the remaining of this section that such an architecture is
sufficient for our purposes. We checked that deeper neural networks do not
provide improvements in the final results, in particular no significant change
can be seen in the autocorrelation times.

###  5.2 Dynamics of the neural network during training

$\beta=0.3$: |  |  |
---|---|---|---
|  |  |
$\beta=0.6$: |  |  |
|  |  |
Figure 2: Dynamics of the neural network training using the KL divergence loss function Eq. (30). For two inverse temperatures: $\beta=0.3$ (upper row) and $\beta=0.6$ (lower row) we provide snapshots after 1, 300, 4000 and 30000 epochs of $\log q_{\theta}$ vs $\log p$. It is evident from the plots that the training is more efficient for smaller values of $\beta$: dots representing states position themselves on $\log q_{\theta}=\log p$ line (green) much faster. The yellow horizontal line shows a uniform probability distribution of $p(\mathbf{s})=2^{-16}$. $\beta=0.3$: |  |  |
---|---|---|---
|  |  |
$\beta=0.6$: |  |  |
|  |  |
Figure 3: Dynamics of the neural network training using the KL divergence loss
function Eq. (30). For two inverse temperatures: $\beta=0.3$ (upper row) and
$\beta=0.6$ (lower row) we provide snapshots after 1, 300, 4000 and 30000
epochs of histograms of $\log w_{i}$.

The small size of the $4\times 4$ Ising model allows us to monitor the
dynamics of the neural network during the training process. All loss functions
have a common, unique minimum where the probabilities provided by the neural
network, $\log q_{\theta}$, are equal to the Boltzmann probabilities, $\log
p$. By plotting the former against the latter at several values of the
training time we can monitor how the differences between the two evolve during
the training. In the presented case we used the KL divergence loss function
Eq. (30). Other loss functions show a similar quantitative behaviour, with
differences at the intermediate stages of the training associated to the
different weighting of contributions from particular configurations. At the
initial moment the weights in the neural network are random and all
probabilities are approximately equal. They form a horizontally elongated
cloud along the yellow line shown in Fig. 2. The yellow line corresponds to a
uniform probability distribution, $q_{\theta}(\mathbf{s})=2^{-16}$. When $\log
q_{\theta}=\log p$ the data points form a straight line, as can be seen in the
last stages of the training shown in the most right hand side plots in Fig. 2.
The two plots in-between are the most interesting ones as they capture the way
the algorithm improves on $q_{\theta}$. As expected, first the configurations
which have the largest probabilities $\log p$ (upper right corner of the
figures) are tuned to the desired values and position themselves on the
straight line. Once they have the correct values the training improves the
approximation of less probable configurations. Please note that at the same
time the approximated probabilities for configurations which are very rare are
also improved, although it is very unlikely that they have appeared in any
batch. This is due to the generalisation properties of the network.

The same information but plotted in somewhat different manner is provided in
Fig. 3 where we show snapshots of the histograms of the log of importance
ratios $\log w$. In the ideal case, the neural network would provide $\log
q_{\theta}$ which are equal to $\log p$ and all importance ratios would be
equal to 1. Indeed, as the training proceeds the histogram is more and more
peaked at $\log 1=0$. It is also evident that it is much easier for the neural
network to learn the probability distribution $p$ at $\beta=0.3$ than at
$\beta=0.6$.

###  5.3 Variational free energy

Figure 4: Free energy $F$ for system of size $4\times 4$ at $\beta=0.6$. We
show value estimated in the batch $\hat{F}_{q}$ (red dots) and value $F_{q}$
calculated using all states (black curve). They are compared with analytical
result for $F$ (green horizontal line). It has to be emphasized that the red
data obtained from the batch averages follows closely the full estimate of the
free energy plotted with the black line.

Now we turn our attention to the variational estimation of the free energy
given by Eq. (29). This is a particular observable for which we search for a
variational approximation from above using the setup we have described in
previous section. Note, that the results presented below are not outcomes from
a Monte Carlo simulation. The discussion below will serve two purposes. On one
hand, we demonstrate how accurately the correct probability distribution can
be reproduced from the point of view of a physically relevant observable. On
the second hand, this benchmark observable can be used to select the simplest
neural network architecture which will be sufficient for our purposes. In
Fig.4 we show the free energy evaluated in three different ways as a function
of the training epoch at $\beta=0.6$. We used a small neural network
architecture – 2 dense layers with autoregressive connections. Each layer
contains $16$ (i.e. number of spins) neurons. We use the $Z_{2}$ symmetry of
the system (see also Section 7.1). Different data points in Fig. 4 correspond
to: Eq. (29) which is the full $F_{q}$ (black curve), the batch estimate
$\hat{F}_{q}$ which is Eq. (29) with replacement Eq. (42) (red dots) and the
value of $F$ calculated analytically (green horizontal line) – see Ref. [13]
or [8] for formulas. The convergence of $F_{q}$ and $\hat{F}_{q}$ to the exact
value $F$ confirms that such a very small neural network architecture was
enough for this system size. Additionally, we can remark that two ways of
calculating the free energy: $F_{q}$ and $\hat{F}_{q}$ agree very well. Note
also that at early epochs $F_{q}$ is smaller than true $F$ because
$\beta$-annealing was used. A similar test performed at other values of
$\beta$ confirmed that the neural network architecture is sufficient. When
considering larger system size we will keep the number of layers equal to 2
and only adjust their width to match the increasing number of spins.

###  5.4 Autocorrelation times

Figure 5: Four versions of (integrated) autocorrelation time (defined in
section 3) for system of size $4\times 4$. Left: autocorrelation times during
the training process for $\beta=0.6$. Right: autocorrelation times as function
of $\beta$, they were measured at epoch 15000, when networks were fully
trained. The system size $4\times 4$ allows us to compare the theoretical
integrated autocorrelation time with the other definitions. The behaviour of
autocorrelation times for a large system size $16\times 16$ is shown in Fig.8
in which case theoretical autocorrelation time $\tau_{L}^{int}$ is not
available.

The next quantity which we can discuss at this small system size are the
different integrated autocorrelation times defined in section 3 by equations
(19), (21), (23) and (26). We start by showing how the autocorrelation time
changes during the training process. The left panel of Fig. 5 shows the
evolution at $\beta=0.6$ as a function of training epochs. The most striking
is the black curve corresponding to the theoretical integrated autocorrelation
time $\tau_{L}^{int}$, Eq. (19). It provides the reference scale for the worst
possible behaviour. It is interesting that after a period of very pronounced
values larger than $10^{4}$ it decreases and stabilizes around 10 when the
training is finished. The transition time coincides with the time at which the
variational free energy converges towards the exact value, see Fig. 4, i.e.
around 8000 epochs. Similar behaviour can be seen for the other estimators,
they all show a correlated drop in this period of training.

The large discrepancy between $\tau_{L}^{int}$ and the other three curves has
two sources. The large value of $\tau_{L}^{int}$ is driven by some
underestimated configuration which has a very small probability to appear in
the batch, therefore, it will not contribute to $\hat{\tau}^{int}_{B}$. The
configurations which contribute to $\hat{\tau}^{int}_{B}$ are those which are
generated in the batch and these are the same configurations which are used to
improve $q_{\theta}$. Hence, $\hat{\tau}^{int}_{B}$ is estimated from
configurations whose importance ratios are closer to one, leading to smaller
autocorrelation time. The other source is related to the fact that
$\hat{\tau}^{H,int}_{F}$ and $\hat{\tau}^{H,int}$ are calculated for a
specific observable, in our case the energy, and hence may couple differently
to different modes of the Markov chain. The discrepancy with $\tau_{L}^{int}$
signals that the contribution of the slowest mode is small in the case of such
physical observables as energy or magnetization. As far as the statistical
uncertainties of such physical observables are concerned, it will be
$\hat{\tau}^{H,int}$ and $\hat{\tau}^{H,int}_{F}$ which matter. It is worth
noticing that $\hat{\tau}^{int}_{B}$ is very close to $\hat{\tau}^{H,int}_{F}$
and $\hat{\tau}^{H,int}$, hence it provides a reasonable estimator of the
autocorrelation time in the Markov chain without making reference to any
specific observable. The increased fluctuations of $\hat{\tau}^{int}_{B}$ are
present because at each epoch new states in the batch are drawn and different
importance ratios appear.

In right panel of Fig. 5 we show the dependence of the various definitions of
the autocorrelation time estimated after 15000 epochs of training as a
function of $\beta$. From the traditional Metropolis algorithm we expect a
sharp maximum around the ordered/disordered phase cross-over inverse
temperature around $\beta\approx 0.4$, however no rise of autocorrelation time
is observed in this region. What we see instead is a strong rise of
$\tau_{L}^{int}$ which signifies that there are configurations for which
$q_{\theta}$ provides an increasingly badly underestimated approximation to
the target probability. Such configurations have a very small probability of
being generated and so they did not appear in any batch hence do not
contribute to any other autocorrelation time, $\hat{\tau}^{H,int}$,
$\hat{\tau}^{H,int}_{F}$ and also $\hat{\tau}^{int}_{B}$ are always of the
order of 1.0.

##  6 Medium lattice: loss functions and critical slowing down

We now move to describe results from larger lattices. We study the Ising model
on a lattice with extent $L=16$. For this system size the eigenvalues given by
the analytic formula derived by Liu, Eq. (12), cannot be evaluated directly.
Instead, one has to rely on quantities which can be estimated as batch
averages. We showed how the two, the exact one and the batch-approximated one,
correlate in a smaller system in the previous section.

We consider different loss functions and study the behaviour of several
observables as a function of training epochs to visualise some aspects of the
neural network training process. Again we keep the simplest neural network
architecture, with two dense layers of width equal to the number of spins on
the lattice: 256. We checked that deeper architectures do not provide an
improvement of the results, but they increase the necessary number of epochs
in the training phase.

###  6.1 Comparison of loss functions

|
---|---
|
|
Figure 6: Learning process for $16\times 16$ system at $\beta=0.6$ using three loss functions. Different observables has been shown as function of training epochs. Functions shown on the plots averages of 16 learning process in order to exclude effects of random net initialization. |
---|---
|
|
Figure 7: The same as Fig. 6 but for $\beta=0.3$

We have mentioned in Section 4 that the different loss functions are sensitive
to different regions of the support of $q_{\theta}$ probability distribution.
In order to try to quantify the effect of these differences on the observables
accessible from the level of a simulation we now take a closer look at six
quantities: the free energy $\hat{F}_{q}$, $\widehat{ESS}$, acceptance ratio
in NMCMC and three definitions of the integrated autocorrelation times. We
show all of them as a function of training epochs in Figs. 6 and 7, with each
of them calculated with the three loss functions introduced earlier. In order
to have a better perspective we compare the situation at $\beta=0.6$ (Fig. 6)
and $\beta=0.3$ (Fig. 7). All data points have been averaged over 16
simulations with different, random initial weights of the neural network. We
also stress that at $\beta=0.3$ we impose the translational symmetry (to be
described and discussed in Section 7.1).

In the first panel of Fig. 6 we present the variational batch estimate of the
free energy. The green line denotes the exact result and one clearly sees that
the Kullback–Leibler divergence loss function is the most effective, i.e.
provides the best approximation at a given number of training epochs. The
$\hat{D}_{ESS}$ loss function fails to converge to the expected value for this
larger system. The authors of Ref. [20] also note in the context of
normalizing flows that this function is unsuitable for training. The same
behaviour of the discussed loss functions is reproduced in other observables.
In other words, the better approximated is the free energy, the larger is the
$\widehat{ESS}$ coefficient, the larger is the acceptance ratio and the
smaller is the integrated autocorrelation time $\hat{\tau}^{H,int}$. Also
other estimates of the autocorrelation time, $\hat{\tau}^{H,int}_{F}$ and
$\hat{\tau}^{int}_{B}$ exhibit the same ordering. Analogous conclusions can be
drawn from Fig. 7 where we present simulations performed at $\beta=0.3$. Note
that the $\hat{D}_{ESS}$ loss functions is unable to show any improvement over
training period.

We note that in all panels of Fig. 6 the observables behave monotonously with
the training time. Some of them show a scatter related to the variations
between consecutive batches, however this does not affect the general trend.
In particular, we do not see any period of stagnation of the training process.

When comparing the two most effective loss functions, $\hat{D}_{KL}$ and
$\hat{D}_{\lambda_{1}}$, we observe that the failure to reproduce the correct
value of the free energy by $\hat{D}_{\lambda_{1}}$ by $0.5$ ‰ corresponds to
an increase of the integrated autocorrelation time $\hat{\tau}^{H,int}$ by a
factor 2.

###  6.2 Critical slowing down

Figure 8: Left: comparison of the autocorrelation time $\hat{\tau}^{H,int}$
between different algorithms for system of size $16\times 16$ as functions of
$\beta$. Right: different autocorrelation times in NMCMC measured at epoch
30000, when network was fully trained.

The central issue of the MCMC simulations is the critical slowing down – the
increase of autocorrelation time when system is simulated closer to the phase
transition. In the traditional approach to the simulations of the Ising model
with algorithms built upon local updates this can be very pronounced as is
shown in the left panel of Fig. 8. For the sake of comparison, we show there
the integrated autocorrelation time $\hat{\tau}^{H,int}$ measured for the
Markov chain obtained using three algorithms: the NMCMC, the original
Metropolis algorithm [1] and the Wolff’s cluster algorithm [11, 12]. The
latter is a well-known case of cluster algorithm where the critical slowing
down was completely eliminated. In order to compare the efficiencies of the
Wolff and Metropolis algorithm we measured the energy after one sweep (a sweep
consists of $L^{2}$ proposals so that each spin of the lattice is updated). In
the case of cluster algorithm one sweep was defined as $L^{2}$ updates of
clusters initiated at each spin of the lattice.

The autocorrelation time $\hat{\tau}^{H,int}$ in the NMCMC algorithm (magenda
curve) clearly exhibits a increase around $\beta\approx 0.44$ with the maximum
value close to 2. Still, such an increase is much less pronounced than the one
of the Metropolis algorithm (red curve). Although, from the presented
comparison the cluster algorithm (black curve) is clearly preferred, one
should remember that there exist relevant physical system, for instance
lattice field theories with fermionic degrees of freedom, for which a cluster
formulation is unknown, whereas there are no limitations in principle in
implementation of the NMCMC approach for other systems.

In the right panel of Fig. 8 we show, for NMCMC, the dependence of various
estimates of the integrated autocorrelation time on $\beta$. A similar result
was discussed in the context of the small system size in Fig. 5 where no
significant dependence on $\beta$ was seen for these autocorrelation times.
Here we clearly observe a maximum around $\beta\approx 0.44$ for all three
quantities.

##  7 Reduction of autocorrelation time

In the NMCMC approach the autocorrelation times can be reduced by improving
training hence reducing the difference between $q_{\theta}$ and $p$. Here we
shall discuss two ideas which help obtaining this goal: i) making use of the
symmetries of the model and ii) reducing effective number of spins using
nearest neighbour interaction property. We shall test those two approaches
using larger system size, i.e. $16\times 16$ spin lattices.

###  7.1 Impact of symmetry on training efficiency and autocorrelation times

In the NMCMC algorithm discussed so far we did not discuss any of the
symmetries of the simulated system. We define symmetries as global
transformations of configurations which do not change their energies. It
immediately follows that configurations related by a symmetry have equal
target probabilities. In order to be more precise, consider a finite group
consisting of $I$ elements
$\\{\mathcal{S}_{0},\mathcal{S}_{1},\ldots,\mathcal{S}_{I-1}\\}$, where
$\mathcal{S}_{0}\equiv 1$ is identity transformation
$\mathcal{S}_{0}\mathbf{s}=\mathbf{s}$. The transformations $\mathcal{S}$
leave the energy unchanged, i.e.
$H(\mathbf{s})=H(\mathcal{S}_{1}\mathbf{s})=\ldots=H(\mathcal{S}_{I-1}\mathbf{s})$.
Symmetries are not present in the conditional probabilities parametrized by
the neural network and hence there is no obvious way of imposing them on the
output of the neural network. The approach suggested in Ref.[7] was to force
the network to treat equivalent configurations as equally probable by
symmetrizing their approximated probabilities within the batch. For that
purpose, for each configuration $\mathbf{s}$ proposed by the neural network we
replace its probability $q(\mathbf{s})$ by the following average:

$q_{\theta}(\mathbf{s})\rightarrow\frac{1}{I}\sum_{i=0}^{I-1}q_{\theta}(\mathcal{S}_{i}\mathbf{s})$
(45)

where $q_{\theta}(\mathcal{S}_{i}\mathbf{s})$ are obtained by applying network
to the transformed spin configuration. This induces an additional numerical
cost because for each configuration from the batch we need to call the neural
network $I-1$ more times to estimate the probabilities of transformed
configurations. However, note that already in order to generate the original
configuration $\mathbf{s}$ we have to evaluate $L^{2}$ conditional
probabilities, meaning that we have called the neural network $L^{2}$ times.
Hence, if $I\ll L^{2}$ the additional cost of imposing the symmetry is
subleading, whereas the benefits may be significant, as we discuss below.

In the case of the Ising model we consider two different types of symmetry:

  * 1.

the $Z_{2}$ symmetry was already discussed in Ref.[7]. It amounts to the
statement that flipping all the spins of a configuration produces a
configuration with the same energy. In this case $I=2$ and

$\mathcal{S}^{Z_{2}}_{1}\mathbf{s}_{x,y}=-\mathbf{s}_{x,y}.$ (46)

Eq.(45) becomes
$q_{\theta}(\mathbf{s})\rightarrow\left[q_{\theta}(\mathbf{s})+q_{\theta}(-\mathbf{s})\right]/2$,

  * 2.

periodic boundary conditions in both directions of the plane allow for
existence of translational symmetries which we denote $T_{x}$ and $T_{y}$,
each composed of $L$ generators, $I=L$, defined as

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

In the following we only implement the $T_{y}$ symmetry. Likewise, Eq.(45) in
this situation becomes

$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,$
(48)

Symmetries ${Z_{2}}$ and ${T_{y}}$ are independent and can be imposed
separately or simultaneously. We show their effects in Figs. 9 and 10
respectively for $\beta=0.6$ and $\beta=0.3$

|
---|---
|
|
Figure 9: Four quantities (free energy, accept ratio in NMCMC, $\widehat{ESS}$ and integrated autocorrelation time $\hat{\tau}^{H,int}$) as functions of epochs during training for system $16\times 16$ and $\beta=0.6$. For each quantity we compare the quality of learning process for different symmetries. Free energy and accept ratio we show zoomed in order to distinguish between the most efficient approaches.  |
---|---
|
|
Figure 10: The same as Fig. 9 but for $\beta=0.3$.

by discussing four observables: free energy $\hat{F}_{q}$, acceptance ratio in
NMCMC, $\hat{E}SS$ and $\hat{\tau}^{H,int}$. Each figure contains five sets of
data: results from the NMCMC simulation without the imposition of any of the
system symmetries, NMCMC simulation with the $Z_{2}$ symmetry alone, NMCMC
simulation with the $T_{y}$ symmetry alone, NMCMC simulation with both $Z_{2}$
and $T_{y}$ symmetries imposed simultaneously and finally a data set obtained
from a NMCMC simulation where the heathbath updates were incorporated on a
chessboard lattice. The details of the latter are described in the next
subsection. Figures showing the free energy contain additionally a green line
at the exact value. All plots show the evolution of the four observables as a
function of the training period, so we expect convergence of the different
data sets at the right edges of the figures.

Fig. 9 shows the results obtained at $\beta=0.6$ when the system is expected
to be in the ordered phase. From the first plot one readily concludes that the
imposition of $Z_{2}$ symmetry is crucial. Data sets corresponding to
simulation without this symmetry converge very slowly or do not converge at
all towards the exact value of the free energy. Only when we include the
$Z_{2}$ symmetry the algorithm is able to reproduce the expected value which
is clearly visible on the zoomed plot. Other observables exhibit a similar
behaviour: acceptance rate is of the order of 40% without the imposition of
$Z_{2}$ symmetry, but once it is included it jumps above 90%.

The situation is different at $\beta=0.3$ which we show in Fig. 10. At this
temperature the system is expected to be in the disordered phase and the role
of $Z_{2}$ and $T_{y}$ symmetries is reversed. In this situation the
imposition of the $Z_{2}$ symmetry does not produce any visible improvement:
for all observables the data set obtained from the simulation without any
symmetry and with the $Z_{2}$ symmetry only are barely distinguishable. It is
now the translational symmetry which improves the performance of the
algorithm. Let us also note that the training is much more efficient than it
was at $\beta=0.6$ as all four observables reach satisfactory values already
after 10000 epochs.

Comparing the results shown in Figs. 9 and 10 one concludes that the
imposition of system symmetries is pivotal in order to make the training
process more efficient. Additionally, we observe that the different symmetries
play different role in the various regions of the phase space of the system,
and hence one should implement as many of them as possible. Let us highlight
again that in the case of two-dimensional systems, as the system size grows,
the number of invocations of the neural network scales like $L^{2}$ whereas
the cost of the symmetries depends on the number of symmetry transformations
that are used to average the approximated probabilities.

###  7.2 Chessboard factorization

Figs. 9 and 10 contained another data set which was not commented in the
previous section and which we explain below. The chessboard factorization
which was used to generated this data can be applied to models with only
nearest-neighbours interactions. One has to divide all spins into two species,
odd and even, according to the parity of the sum of their coordinates on the
lattice, this results in a ”chessboard” pattern. From the locality of the
interaction it turns out that the total probability of that spin from an even
site depends exclusively on the state of the four neighbouring spins from the
odd sites. In order words, if all spins at odd sites have been fixed, the
spins at even sites can be generated from the exact, local Boltzmann
probability distribution.

In the present context the above formulation means that one can reduce the
number of spins which have to be set by the neural network by a half. Not only
this reduces significantly the number of invocations of the neural network,
but also the size of the neural network is reduced and more importantly half
of the spins can be generated independently from each other from the target
probability distribution further decorrelating the system.

The formula Eq. (4) is now modified:

$\displaystyle
q_{\theta}(\mathbf{s})=\prod_{i=1}^{N/2}q_{O}(s^{i}|s^{i-1}\dots
s^{1})\prod_{i=N/2+1}^{N}q_{E}(s^{i},n(s^{i})),$ (49)

where $q_{O}(s^{i}|s^{i-1}\dots s^{1})$ is a probability of $i$-th _odd_ spin
to be $s^{i}$ when other _odd_ spins $s^{i-1}\dots s^{1}$ are fixed and
$q_{E}(s^{i},n_{A}(s^{i})$ is probability of $i$-th _even_ spin to be $s^{i}$
when it’s neighbours $n(s_{i})$ are fixed (note that even spins have only odd
neighbours). The procedure to get the full spin state is the following: i) run
forward network to get odd spins, ii) using values of fixed odd spins
calculate probabilities $q_{E}$ of all even spins, iii) draw the even spins
using calculated $q_{E}$. Probability of $i$-th even spin to be $+1$ or $-1$
is equal:

$\displaystyle q_{E}(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},$ (50)

where the sum is performed over nearest neighbours $n(s^{i})$ of spin $s_{i}$.

As the neural network is used to generate only half of spins in full state we
expect that training is faster and more effective when chessboard
factorization is used. In Figures 9 and 10 we show (light blue diamonds)
observables during the training when this factorization is applied, no
symmetries involved here (so should be compared with magenta triangles). We
observe dramatic difference in performance of the network - chessboard
factorization improves significantly all measures of training quality,
comparable to that of the symmetries imposition.

##  8 Summary and outlook

In this work we have investigated, from many perspectives, the neural Markov
Chain Monte Carlo algorithm for simulating statistical systems with a finite
number of discrete degrees of freedom. We quantified its performance by
measuring various estimators of integrated autocorrelation times and discussed
their behaviour as a function of inverse temperature $\beta$ and neural
network training time. We have verified that the NMCMC algorithm is mildly
affected by the critical slowing down close to the phase transition, i.e. the
integrated autocorrelation time does not increase considerably in this region.
Of course, this is only true if the neural network was trained well enough so
that the encoded probability distribution is sufficiently close to the target
one. We suggested some measures of the quality of the training based on
importance ratios. The latter are the building blocks of the analytic
expressions for the eigenvalues of the Markov chain transition matrix and
hence are used to describe the autocorrelations in the Metropolized
Independent Sampler class of algorithms. Comparing with the local Metropolis
algorithm we stress that the critical slowing down of the Metropolis algorithm
can be explained by the growing size of spin clusters in the parameter region
around the phase transition, hence has a physical interpretation. In the case
of the NMCMC algorithm the non-zero autocorrelation times are due to an
imperfect approximation of the Boltzmann distribution by $q_{\theta}$.
Therefore the phenomenon of the critical slowing down, if any, stems from the
dynamics of the neural network training. Furthermore the configurations with
very small $p$ are the ones for which the network will have biggest difficulty
for getting $q_{\theta}$ right. This can lead to very large relative errors
and hence very large values of $\tau_{L}$. As those configurations have
$q_{\theta}(\mathbf{s})\ll p(\mathbf{s})$ they will be proposed very rarely if
ever. However if they do appear they will repeated approximately
$w(\mathbf{s})$ times. At this moment it is unclear to us if this may have an
impact on the results of simulations, as we have not encountered such
situation in practice.

As far as the loss functions are concerned we have constructed a loss function
directly proportional to the largest eigenvalue of the transition matrix of
the Markov chain. Although its performance is inferior to the previously
studied $\hat{D}_{KL}$ loss function, it provides an interesting example of a
loss function which behaves differently when the neural network training is
implemented with the REINFORCE algorithm or with the gradient learning. This
is in contrast to $\hat{D}_{KL}$ loss function that is the same in both
approaches. Hence, in further studies $\hat{D}_{\lambda_{1}}$ may be used to
investigate different approaches to neural network training. As a third
alternative, we checked the $\hat{D}_{ESS}$ loss function discussed in Ref.
[20], and confirmed that its efficiency is very poor in the studied context.

We have investigated the impact of symmetries on the training process. It was
already noted in Ref. [7] that the $Z_{2}$ symmetry is crucial in order to
correctly reproduce the free energy of the system. In our detailed study we
observed that this is true only in a specific range of inverse temperatures
$\beta$. For instance, for small $\beta$ the $Z_{2}$ symmetry does not improve
the quality of training. Instead, another symmetry turns out to be more
important, namely the translational symmetry, which in turn was not as
efficient at large values of $\beta$ as the $Z_{2}$ symmetry. Hence, we
conclude that in order to efficiently train the neural network in the full
range of $\beta$ it may be advisable to impose both, or all accessible,
symmetries simultaneously.

Our third observation concerns the chessboard factorization, which turns out
to considerably reduce both the width of the neural network and the training
time. We find that its effects are as important as the effects of the
symmetries and therefore it may be advantageous to incorporate it in the
algorithm in the full range of $\beta$.

In summary we note that the improvements proposed and discussed in this work
also apply to the variational estimation of free energy, as proposed in Ref.
[7].

We leave as an open question and as a next step in the investigation of NMCMC
algorithms, how many of the above improvements can be applied in the context
of systems with continuous degrees of freedom.

## Acknowledgement

Computer time allocation ’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) Grant No. 2019/32/C/ST2/00202.

## Author contributions

All Authors have equally contributed to all stages of this work and of the
manuscript preparation.

## Declaration of competing interest

The Authors declare no competing interests.

## References

  * [1] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953.
  * [2] K. Binder and D. Heermann, Monte Carlo Simulation in Statistical Physics: An Introduction.  Springer, 2019.
  * [3] A. Laio and F. L. Gervasio, “Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science,” Reports on Progress in Physics, vol. 71, p. 126601, Dec. 2008\.
  * [4] G. Bussi and A. Laio, “Using metadynamics to explore complex free-energy landscapes,” Nature Reviews Physics, vol. 2, p. 200, 2020.
  * [5] F. Fucito, C. Rebbi, and S. Solomon, “Finite Temperature QCD in the Presence of Dynamical Quarks,” Nucl. Phys. B, vol. 248, pp. 615–628, 1984\.
  * [6] M. G. Endres, R. C. Brower, W. Detmold, K. Orginos, and A. V. Pochinsky, “Multiscale monte carlo equilibration: Pure yang-mills theory,” Phys. Rev. D, vol. 92, p. 114516, Dec 2015.
  * [7] D. Wu, L. Wang, and P. Zhang, “Solving Statistical Mechanics Using Variational Autoregressive Networks,” Phys. Rev. Lett., vol. 122, p. 080602, Mar. 2019.
  * [8] K. A. Nicoli, S. Nakajima, N. Strodthoff, W. Samek, K.-R. Müller, and P. Kessel, “Asymptotically unbiased estimation of physical observables with neural samplers,” Phys. Rev. E, vol. 101, p. 023304, Feb. 2020.
  * [9] B. McNaughton, M. V. Milošević, A. Perali, and S. Pilati, “Boosting Monte Carlo simulations of spin glasses using autoregressive neural networks,” Phys. Rev. E, vol. 101, p. 053312, May 2020.
  * [10] D. Wu, R. Rossi, and G. Carleo, “Unbiased Monte Carlo Cluster Updates with Autoregressive Neural Networks,” arXiv e-prints, p. arXiv:2105.05650, May 2021.
  * [11] U. Wolff, “Collective monte carlo updating for spin systems,” Phys. Rev. Lett., vol. 62, pp. 361–364, Jan 1989.
  * [12] U. Wolff, “Comparison between cluster monte carlo algorithms in the ising model,” Physics Letters B, vol. 228, no. 3, pp. 379–382, 1989.
  * [13] A. E. Ferdinand and M. E. Fisher, “Bounded and inhomogeneous ising models. i. specific-heat anomaly of a finite lattice,” Phys. Rev., vol. 185, pp. 832–846, Sep 1969.
  * [14] J. S. Liu, “Metropolized independent sampling with comparisons to rejection sampling and importance sampling,” Statistics and Computing, vol. 6, no. 2, pp. 113–119, 1996.
  * [15] W. K. Hastings, “Monte Carlo Sampling Methods Using Markov Chains and Their Applications,” Biometrika, vol. 57, pp. 97–109, 1970.
  * [16] L. Tierney, “Markov Chains for Exploring Posterior Distributions,” The Annals of Statistics, vol. 22, no. 4, pp. 1701 – 1728, 1994.
  * [17] G. Roberts and J. Rosenthal, “Markov-Chain Monte Carlo: Some Practical Implications of Theoretical Results,” The Canadian Journal of Statistics / La Revue Canadienne De Statistique, vol. 26, pp. 5–20, 1998.
  * [18] A. Sokal, “Monte carlo methods in statistical mechanics: Foundations and new algorithms,” in Functional Integration: Basics and Applications (C. DeWitt-Morette, P. Cartier, and A. Folacci, eds.), (Boston, MA), pp. 131–192, Springer US, 1997.
  * [19] A. B. R.S. Sutton, Reinforcement Learning: An Introduction. Cambridge, MA: MIT Press, 2018.
  * [20] M. S. Albergo, D. Boyda, D. C. Hackett, G. Kanwar, K. Cranmer, S. Racanière, D. Jimenez Rezende, and P. E. Shanahan, “Introduction to Normalizing Flows for Lattice Field Theory,” arXiv e-prints, p. arXiv:2101.08176, Jan. 2021.


# Experimental method for perching
flapping-wing aerial robots
††thanks:  The research is funded by the European Project GRIFFIN ERC Advanced
Grant 2017, Action 788247.

Raphael Zufferey  GRVC Robotics Lab,
University of Seville
Seville, Spain, and
EPFL, Switzerland
<EMAIL_ADDRESS>
Daniel Feliu-Talegón  GRVC Robotics Lab
University of Seville
Seville, Spain
<EMAIL_ADDRESS>
Saeed Rafee Nekoo  GRVC Robotics Lab
University of Seville
Seville, Spain
<EMAIL_ADDRESS>
Jose-Angel Acosta  GRVC Robotics Lab
University of Seville
Seville, Spain
<EMAIL_ADDRESS>
Anibal Ollero  GRVC Robotics Lab
University of Seville
Seville, Spain
<EMAIL_ADDRESS>

###### Abstract

In this work, we present an experimental setup and guide to enable the
perching of large flapping-wing robots. The combination of forward flight,
limited payload, and flight oscillations imposes challenging conditions for
localized perching. The described method details the different operations that
are concurrently performed within the 4 second perching flight. We validate
this experiment with a 700 g ornithopter and demonstrate the first autonomous
perching flight of a flapping-wing robot on a branch. This work paves the way
towards the application of flapping-wing robots for long-range missions, bird
observation, manipulation, and outdoor flight.

##  I Introduction

Figure 1: Side view of a typical perching flight experiment. The robotic birds
are accelerated by a launching system to 4 m/s. The birds then separate from
the launcher and start the controlled flight phase, where the altitude, pitch,
yaw and speed is regulated by an autopilot through elevator, rudder and
flapping frequency actions. A perching robot is shown in the following
positions: A. At the end of the launcher B. In controlled flapping flight C.
Shortly before contact, with the misalignment compensation system active and
D. Perched with its bistable claw, after the flight.

Flapping-wing flight offers significant benefits over rotor-based propulsion
ranging from noise reduction, increased safety, and potentially higher
maneuverability and efficiencies [1]. While these traits could grant safe
flight around structures, significant challenges remain to land large
flapping-wing robots on perches and branches.

Indeed, the combination of forward flight requirement and vertical center of
mass oscillations arising in large flapping-wing robots vastly increases the
difficulty to perch on objects. Achieving this task - as we often see birds do
seamlessly everyday - is complex in terms of the control required in stall
conditions and in terms of grasping appendage capability. In a robot, this
would require a vehicle capable of carrying a claw system that can hold the
robot upright for manipulation, a vision system to identify the branch and a
control algorithm that can reliably bring the robot there with minimal
remaining velocity. We propose to solve the perching maneuver in a controlled
context as a significant step towards vision-based outdoor perching. Here, we
would like to show a flight experiment to demonstrate autonomous perching of
flapping robots in controlled conditions. Complete information about the
design, modeling and resulting flight performance of the robotic bird can be
found in [2].

##  II Experiment

Demonstration and validation of this experimental method is performed with a
large 700 g, 150 cm-wingspan ornithopter whose design stems from a previous
robot class described in [3]. Successful perching flights were achieved
through a series of state phases, visually represented in Fig. 1 and described
in this work.

###  II-A Reaching flight velocity

Large ornithopters are not capable of hovering in windless conditions, due to
unfavorable scaling. Such robots can only sustain flight with forward velocity
and therefore require a process to reach this velocity initially. While one
can certainly achieve this velocity via a hand throw, this results in a wide
range of initial conditions. We instead propose an automated launcher system,
which boasts sufficient precision for repeatable flights in space-constraints,
indoor conditions, see Fig. 1.A.

The system features a 2 m long double rail which can accelerate a robot to 10
m/s. In our validation scenario, speeds of 4 m/s are sufficient to initiate
flight without losing altitude and without excessive speed. The robotic birds
are held at a perpendicular offset of 50 cm, sufficient to avoid any
collisions between the wing, tail and structure. Launcher mounting of robot
occurs at the trailing edge, to remain as close as possible to the center of
mass. A 20 degree angle-of-attack at launch is selected. Results show that
launches are consistent, both in direction and speed, allowing us to
repeatedly perform controlled flights to a point location at a 14 m distance
from the launcher.

###  II-B Controlled flapping flight

At the end of the rail, the robot’s position and attitude are captured by the
motion capture system and start being wirelessly streamed to the robot’s
flight companion computer. The ornithopter enters a controlled flapping flight
phase, see Fig. 1.B. This is handled by the onboard triple loop pitch-yaw-
altitude controller. This flight control system needs to handle stable flight
with the appendage payload, misalignment compensation system and all the
onboard electronics required for flight and perch.

The forward flight velocity is a key parameter for perching. It should be kept
low to minimize impact forces and increase reaction time before impact. To
achieve this, the pitch angle of the robot needs to be high. Experiments show
that pitch angles above $40\text{\,}\mathrm{\SIUnitSymbolDegree}$ result in
insufficient speed and consequently loss of lateral control. As such, a pitch
angle of $30\text{\,}\mathrm{\SIUnitSymbolDegree}$is appropriate for perching
flights and is sufficient to fly stably. This set-point is maintained by a
proportional-integral (PI) controller. Consequently, when the pitch reaches
$30\text{\,}\mathrm{\SIUnitSymbolDegree}$, the velocity drops to 2.5-3 m/s,
kept until perching.

The altitude target value is achieved through regulated flapping frequency
control action. Using this method, all flight trajectories reached a 2 m set-
point within 8-12 m of flight distance. The flight controller performs
adequately for different branch altitudes, with a maximum average vertical
error of 16 cm and a lateral error of rarely exceeding 60 cm.

###  II-C Branch approach

The claw dimension is smaller than the flight position accuracy. As this is
insufficient to reliably touch the target, we propose an correction method.
The leg of the robotic bird is articulated at the elbow level, permitting
semi-vertical motion of the claws. When the robot gets to within 150 cm of the
branch, an optical detection system is enabled, see Fig. 1.C. The line
detection system located on the claw feeds the relative position of the branch
to the leg microcontroller. The leg servo then corrects the angle of the leg
to align it with the branch, at 50 Hz. During this approach phase, the
flapping motion is maintained and the flight controller remains active. The
perching target, a 80 cm-black-painted branch is held at 2 m height at its
center by a $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ aluminum profile which
reaches out through the safety net.

###  II-D Perching

At 20 cm from the branch, the flapping is stopped. At this point, the forward
velocity of the robot is situated between 2.5 and 3 m/s. Grasping is made
possible thanks to a pair of bistable claws which pivot from the open to the
closed state within 25 ms and with a 2 Nm torque. A combination of soft
silicone pads and spikes offers friction forces sufficient to compensate the
torque resulting from center of mass offsets.

Experiments have demonstrated that we can reach the branch location
repeatedly, within the tolerance range of the leg-claw system. The low forward
flight velocity shortly before landing ensures that forces remain below 150 N,
significantly reducing the likelihood of damage. We report that branch
perching was achieved in 6 out of 9 perching flight tests, and that no damage
occurred in any of the successful flights. The results were validated with a
second robot which perched without re-tuning.

## References

  * [1] G. de Croon, “Flapping wing drones show off their skills,” Science Robotics, vol. 5, no. 44, p. eabd0233, 2020.
  * [2] R. Zufferey, J. T. Barbero, D. F. Talegon, S. R. Nekoo, J. A. Acosta, and A. Ollero, “How ornithopters can perch autonomously on a branch,” Nature Communications 13 (1), 7713, 2022.
  * [3] R. Zufferey, J. Tormo-Barbero, M. Mar Guzman, F. J. Maldonado, E. Sanchez-Laulhe, P. Grau, M. Perez, J. A. Acosta, and A. Ollero, “Design of the High-Payload Flapping Wing Robot E-Flap,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3097–3104, 2021.