#  Stochastic optimization: Glauber dynamics versus stochastic cellular
automata

Bruno Hideki Fukushima-Kimura1, Yoshinori Kamijima2, Kazushi Kawamura3 and
Akira Sakai1
1Faculty of Science, Hokkaido University
Kita 10, Nishi 8, Kita-ku, Sapporo, Hokkaido 060-0810, Japan
2 National Center for Theoretical Sciences (NCTS)
No. 1, Sec. 4, Roosevelt Rd., Taipei City 106, Taiwan
3Institute of Innovative Research, Tokyo Institute of Technology
Nagatsuta-cho 4259-J3-30, Midori-ku, Yokohama, Kanagawa 226-8503, Japan
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The topic we address in this paper concerns the minimization of a Hamiltonian
function for an Ising model through the application of simulated annealing
algorithms based on (single-site) Glauber dynamics and stochastic cellular
automata (SCA). Some rigorous results are presented in order to justify the
application of simulated annealing for a particular kind of SCA. After that,
we compare the SCA algorithm and its variation, namely the $\varepsilon$-SCA
algorithm, studied in this paper with the Glauber dynamics by analyzing their
accuracy in obtaining optimal solutions for the max-cut problem on Erdös-Rényi
random graphs, the traveling salesman problem (TSP), and the minimization of
Gaussian and Bernoulli spin glass Hamiltonians. We observed that the SCA
performed better than the Glauber dynamics in some special cases, while the
$\varepsilon$-SCA showed the highest performance in all scenarios.

##  1 Introduction

In several problems involving complex networks, there has been a growing
demand to provide algorithms that can offer optimal solutions for large-scale
combinatorial problems within a relatively short amount of time. Due to the
fact that no polynomial-time algorithms for NP-hard problems are known [8], an
alternative approach is necessary. It follows from the fact that a
considerable amount of problems of practical interest can be mapped into the
problem of finding a ground state of an Ising model (e.g., [2, 4, 16, 19]),
there are several algorithms being developed recently aiming at approaching
such ground states, see [1, 10, 11, 17, 21].

The robustness of the statistical mechanical framework has received a lot of
attention in recent decades, and a stochastic approach to combinatorial
problems has been widely employed. The approach adopted in this paper111The
material of this paper was partially presented at the 53rd Annual Conference
of the Institute of Systems, Control and Information Engineers (SCI’21) which
was held in October, 2021 continues the one discussed in [7], and we extend
some of its theoretical results and provide further examples. As in the
previous work, our main problem of interest is the problem of determining one
of the ground states (i.e., the points of global minima) of a Hamiltonian
function $H$ for an Ising model on a finite simple graph $G=(V,E)$. Given a
collection of spin-spin coupling constants $(J_{x,y})_{x,y\in V}$ such that
$J_{x,y}=J_{y,x}$, and $J_{x,y}=0$ if $\\{x,y\\}\notin E$, and local external
fields $(h_{x})_{x\in V}$, let us consider the Ising Hamiltonian defined by

$\displaystyle H(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{x,y\in
V}J_{x,y}\sigma_{x}\sigma_{y}-\sum_{x\in V}h_{x}\sigma_{x}$ (1.1)

for each spin configuration $\boldsymbol{\sigma}=(\sigma_{x})_{x\in
V}\in\\{-1,+1\\}^{V}$.

Approximating an Ising Hamiltonian’s ground states can be achieved by applying
simulated annealing based on a Markov chain Monte Carlo spin-flip dynamics.
The most usual example consists of considering a time-inhomogeneous Markov
chain based on a single spin-flip dynamics (such as the Glauber or Metropolis
dynamics [9]) while slowly decreasing its temperature towards zero. Some
theoretical results (e.g., [3, 12]) show that if the temperature is decreased
in time $t$ at a certain speed rate, such a procedure yields an approximation
to the minimum value of $H$. In [17, 21], the authors considered simulated
annealing based on parallel spin-flip dynamics that, differently from the
usual method, allows the system to update multiple spins independently, and
simulations showed that it manages to converge to the ground states
significantly faster compared to the single spin-flip methods. In addition to
these observations, in [6] the authors investigated from the mathematical
point of view a specific kind of stochastic cellular automata (SCA) which was
derived from [18, 20], in order to determine sufficient conditions under which
the convergence to the ground states is guaranteed.

As introduced in [6, 13, 18, 20], the SCA transition kernel at inverse
temperature $\beta\geq 0$ for a collection of (non-negative) pinning
parameters $\boldsymbol{q}=(q_{x})_{x\in V}$ is given by

$\displaystyle P^{\scriptscriptstyle\rm
SCA}_{\beta,\boldsymbol{q}}(\boldsymbol{\sigma},\boldsymbol{\tau})=\prod_{x\in
V}\frac{e^{\frac{\beta}{2}(\tilde{h}_{x}(\boldsymbol{\sigma})+q_{x}\sigma_{x})\tau_{x}}}{2\cosh(\frac{\beta}{2}(\tilde{h}_{x}(\boldsymbol{\sigma})+q_{x}\sigma_{x}))}$
(1.2)

for any spin configurations $\boldsymbol{\sigma},\boldsymbol{\tau}$ in
$\\{-1,+1\\}^{V}$, where the cavity fields $\tilde{h}_{x}$ are defined by
letting $\tilde{h}_{x}(\boldsymbol{\sigma})=\sum_{y\in
V}J_{x,y}\sigma_{y}+h_{x}$. We immediately notice that the right-hand side of
the equation above is expressed in terms of a product of local probabilities
$p_{x,q}(\cdot|\boldsymbol{\sigma})$ given by

$p_{x,q}(s|\boldsymbol{\sigma})=\frac{e^{\frac{\beta}{2}(\tilde{h}_{x}(\boldsymbol{\sigma})+q\sigma_{x})s}}{2\cosh(\frac{\beta}{2}(\tilde{h}_{x}(\boldsymbol{\sigma})+q\sigma_{x}))}$
(1.3)

for $s\in\\{-1,+1\\}$. Therefore, whenever we update the state of the system,
all spins are updated simultaneously and independently according to a
probabilistic rule for each vertex that depends on the current spin
configuration. Such characteristics imply that the system can potentially
transition from any spin configuration to another in a single step,
contrasting with the standard Glauber and Metropolis dynamics. Note that due
to the presence of the pinning parameters in the local transition
probabilities, the SCA tends to flip a smaller amount of spins from certain
configurations, especially when such parameters are taken sufficiently large
and the temperature approaches zero. It follows that, for SCA-based simulated
annealing algorithms, such an effect typically slows down the dynamics, which
may prevent the system from reaching lower energetic configurations, however,
as we observe in this paper (see also [7]), such algorithms still performs
better than Glauber dynamics with regards to the success rate of reaching a
ground state in certain scenarios.

In addition to the SCA defined in (1.2), let us also consider the so-called
$\varepsilon$-SCA. First introduced in [6, 7], the $\varepsilon$-SCA is also a
particular kind of probabilistic cellular automata, but with the advantage of
having attenuated temperature-dependent pinning effects, which encourages a
larger number of spin-flips at low temperatures compared to the SCA. Given the
inverse temperature $\beta\geq 0$ and a parameter $\varepsilon\in(0,1]$, the
transition kernel of the $\varepsilon$-SCA is given by

$P_{\beta,\varepsilon}(\boldsymbol{\sigma},\boldsymbol{\tau})=\prod_{x:\,\tau_{x}=-\sigma_{x}}\big{(}\varepsilon
p_{x}(\boldsymbol{\sigma})\big{)}\prod_{y:\,\tau_{y}=\sigma_{y}}\big{(}1-\varepsilon
p_{y}(\boldsymbol{\sigma})\big{)}$ (1.4)

for every pair $\boldsymbol{\sigma},\boldsymbol{\tau}$ of configurations in
$\\{\\-1,+1\\}^{V}$, where $p_{x}(\boldsymbol{\sigma})$ is defined by

$p_{x}(\boldsymbol{\sigma})=\frac{e^{-\frac{\beta}{2}\tilde{h}_{x}(\boldsymbol{\sigma})\sigma_{x}}}{2\cosh(\frac{\beta}{2}\tilde{h}_{x}(\boldsymbol{\sigma}))}.$
(1.5)

Note that $p_{x}(\boldsymbol{\sigma})$ coincides with
$p_{x,0}(-\sigma_{x}|\boldsymbol{\sigma})$, so this algorithm can be
interpreted in the following way: the spins eligible to be flipped are
randomly assigned with probability $\varepsilon$ independently of each other,
then, only those spins which were assigned will be updated simultaneously and
independently, where the probability of flipping the spin $\sigma_{x}$ is
$p_{x}(\boldsymbol{\sigma})$, while the unassigned spins remain unchanged.
Note that $P_{\beta,\varepsilon}$ can also be written as a product of local
transition probabilities just like in equation (1.2), where each of the
corresponding local probabilities
$\tilde{p}_{x,\varepsilon}(\cdot|\boldsymbol{\sigma})$ is given by

$\tilde{p}_{x,\varepsilon}(s|\boldsymbol{\sigma})=(1-\varepsilon)\delta_{\sigma_{x},s}+\varepsilon
p_{x,0}(s|\boldsymbol{\sigma})$ (1.6)

for each $s\in\\{-1,+1\\}$.

Differently from the SCA and some of the well-established single spin-flip
dynamics, there is still no solid mathematical theory that describes in which
circumstances a simulated annealing algorithm based on the $\varepsilon$-SCA
can be successfully applied such as in Theorem 2.3. Even though, in recent
investigations, e.g., [6, 7], the $\varepsilon$-SCA outperformed the tested
algorithms in all scenarios, managing to reach lower energy configurations and
the ground states at a higher success rate, provided the parameter
$\varepsilon$ was chosen appropriately. In this paper, we provide a natural
extension for the preliminary results from [7], where we derive a new result
concerning an upper bound for the mixing time of the $\varepsilon$-SCA, and
test the algorithms by applying them to a larger class of problems, including
the traveling salesman problem (TSP).

##  2 Mathematical results

In this section, we summarize some of the rigorous results which are known for
the SCA and the $\varepsilon$-SCA. Furthermore, we also provide an upper bound
for the mixing time of the $\varepsilon$-SCA. For the readers interested in
the proofs concerning to the SCA, refer to [6]. The first mathematical
results, Theorems 2.1 and 2.2, show that the mixing time $t_{mix}(\delta)$
(see definition in [15]), that is, the time the system takes to be close to
equilibrium, is bounded above by a quantity which is proportional to $\log|V|$
as long as the temperature is taken sufficiently high. Such results show that
both the SCA and the $\varepsilon$-SCA reach the equilibrium faster than
Glauber dynamics which is known to be at least proportional to $|V|\log|V|$,
see [5, 14].

######  Theorem 2.1 (Mixing time [6]).

For any non-negative pinning parameters $\boldsymbol{q}=(q_{x})_{x\in V}$, if
$\beta$ is sufficiently small so that

$\displaystyle r\equiv\max_{x\in V}\bigg{(}\tanh\frac{\beta
q_{x}}{2}+\sum_{y\in V}\tanh\frac{\beta|J_{x,y}|}{2}\bigg{)}<1,$ (2.1)

then, the mixing time $t_{\rm mix}(\delta)$ satisfies

$t_{\rm
mix}(\delta)\leq\bigg{\lceil}\frac{\log|V|-\log\delta}{\log(1/r)}\bigg{\rceil}.$
(2.2)

######  Theorem 2.2 (Mixing time of the $\varepsilon$-SCA).

For any positive $\varepsilon\in(0,1]$, if $\beta$ is sufficiently small such
that

$\displaystyle r\equiv(1-\varepsilon)+\varepsilon\max_{x\in
V}\bigg{(}\sum_{y\in V}\tanh\frac{\beta|J_{x,y}|}{2}\bigg{)}<1,$ (2.3)

then, the mixing time $t_{\rm mix}(\delta)$ satisfies

$t_{\rm
mix}(\delta)\leq\bigg{\lceil}\frac{\log|V|-\log\delta}{\log(1/r)}\bigg{\rceil}.$
(2.4)

###### Proof.

Following the same strategy as adopted in the proof of [6, Proposition 2.2],
it is sufficient to prove that $\rho_{\scriptscriptstyle\rm
TM}(P_{\beta,\varepsilon}(\boldsymbol{\sigma},\cdot),P_{\beta,\varepsilon}(\boldsymbol{\tau},\cdot))\leq
r$ holds for any configurations $\boldsymbol{\sigma},\boldsymbol{\tau}$
satisfying $|D_{\boldsymbol{\sigma},\boldsymbol{\tau}}|=1$, where
$\rho_{\scriptscriptstyle\rm TM}$ stands for the transportation metric and
$D_{\boldsymbol{\sigma},\boldsymbol{\tau}}=\\{x\in
V:\sigma_{x}\neq\tau_{x}\\}$.

Let us denote the local probabilities
$\tilde{p}_{y,\epsilon}(+1|\boldsymbol{\sigma})$ from equation (1.6) simply by
$p(\boldsymbol{\sigma},y)$, and let us assume that
$D_{\boldsymbol{\sigma},\boldsymbol{\tau}}=\\{x\\}$, that is,
$\boldsymbol{\tau}$ is the configuration
$\boldsymbol{\tau}=\boldsymbol{\sigma}^{x}$ which coincides with
$\boldsymbol{\sigma}$ apart from its value at $x$. Note that
$p(\boldsymbol{\sigma},y)\neq p(\boldsymbol{\sigma}^{x},y)$ only if $y=x$ or
$y\in N_{x}\equiv\\{v\in V:J_{x,v}\neq 0\\}$. Using this as a threshold
function for i.i.d. uniform random variables $(U_{y})_{y\in V}$ on $[0,1]$, we
define the coupling $(X,Y)$ of
$P_{\beta,\varepsilon}(\boldsymbol{\sigma},\cdot)$ and
$P_{\beta,\varepsilon}(\boldsymbol{\sigma}^{x},\cdot)$ as

$\displaystyle X_{y}=\begin{cases}+1&[U_{y}\leq p(\boldsymbol{\sigma},y)],\\\
-1&[U_{y}>p(\boldsymbol{\sigma},y)],\end{cases}$ (2.5)

and

$\displaystyle Y_{y}=\begin{cases}+1&[U_{y}\leq
p(\boldsymbol{\sigma}^{x},y)],\\\
-1&[U_{y}>p(\boldsymbol{\sigma}^{x},y)].\end{cases}$ (2.6)

If we denote the probability measure of this coupling by
$\textbf{P}_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{x}}$ and its expectation
by $\textbf{E}_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{x}}$, then, it
follows that

$\displaystyle\textbf{E}_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{x}}[|D_{X,Y}|]=$
$\displaystyle|p(\boldsymbol{\sigma},x)-p(\boldsymbol{\sigma}^{x},x)|$
$\displaystyle+\sum_{y\in
N_{x}}|p(\boldsymbol{\sigma},y)-p(\boldsymbol{\sigma}^{x},y)|,$ (2.7)

where,

$\displaystyle|p(\boldsymbol{\sigma},x)-p(\boldsymbol{\sigma}^{x},x)|=1-\varepsilon,$
(2.8)

and for $y\in N_{x}$ we have

$\displaystyle|p(\boldsymbol{\sigma},y)-p$
$\displaystyle(\boldsymbol{\sigma}^{x},y)|$
$\displaystyle\leq\frac{\varepsilon}{2}\bigg{|}$
$\displaystyle\tanh\bigg{(}\frac{\beta(\sum_{v\neq
x}J_{v,y}\sigma_{v}+h_{y})}{2}+\frac{\beta J_{x,y}}{2}\bigg{)}$
$\displaystyle-\tanh\bigg{(}\frac{\beta(\sum_{v\neq
x}J_{v,y}\sigma_{v}+h_{y})}{2}-\frac{\beta J_{x,y}}{2}\bigg{)}\bigg{|}.$ (2.9)

Since $|\tanh(a+b)-\tanh(a-b)|\leq 2\tanh|b|$ for any $a,b$, then, we conclude
that

$\displaystyle\rho_{\scriptscriptstyle\rm TM}\Big{(}P^{\scriptscriptstyle\rm
SCA}_{\beta,\boldsymbol{q}}(\boldsymbol{\sigma},\cdot),P^{\scriptscriptstyle\rm
SCA}_{\beta,\boldsymbol{q}}(\boldsymbol{\sigma}^{x},\cdot)\Big{)}$
$\displaystyle\leq\textbf{E}_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{x}}[|D_{X,Y}|]$
$\displaystyle\leq(1-\varepsilon)$ $\displaystyle+\varepsilon\sum_{y\in
N_{x}}\tanh\frac{\beta|J_{x,y}|}{2}\leq r,$ (2.10)

as required. ∎

The next result shows that if we consider the inhomogeneous Markov chain where
the system is updated at time $t$ according to the SCA transition kernel at
inverse temperature $\beta_{t}$ proportional to $\log t$, then, the algorithm
will converge to the uniform distribution concentrated over the ground states
of the Hamiltonian (1.1).

######  Theorem 2.3 (Simulated annealing [6]).

Let us assume that each $q_{x}\geq\lambda/2$, where $\lambda$ is the largest
eigenvalue of the matrix $[-J_{x,y}]_{x,y\in V}$. If we choose
$\\{\beta_{t}\\}_{t\in{\mathbb{N}}}$ as

$\beta_{t}=\frac{\log t}{\Gamma},$ (2.11)

where $\Gamma=\sum_{x\in V}\Gamma_{x}$ and
$\Gamma_{x}=q_{x}+|h_{x}|+\sum_{y\in V}|J_{x,y}|$, then, for any initial
$j\in{\mathbb{N}}$, we have

$\displaystyle\lim_{t\to\infty}\sup_{\mu}\big{\|}\mu P^{\scriptscriptstyle\rm
SCA}_{\beta_{j},\boldsymbol{q}}P^{\scriptscriptstyle\rm
SCA}_{\beta_{j+1},\boldsymbol{q}}\cdots P^{\scriptscriptstyle\rm
SCA}_{\beta_{t},\boldsymbol{q}}-\pi^{\scriptscriptstyle\rm
G}_{\infty}\big{\|}_{\scriptscriptstyle\rm TV}=0.$ (2.12)

It is important to point out that the condition on the values of the pinning
parameters $\boldsymbol{q}=(q_{x})_{x\in V}$ to be at least $\lambda/2$ is
just a sufficient condition related to a more general one originated in [17].
Therefore, such bound still has some room for improvement and the algorithm
may still converge even if we consider smaller values for such parameters. The
importance of letting the parameters $q_{x}$ assume relatively large values
comes from the fact that it prevents the system from flipping many spins at
the same time, which helps to guarantee the convergence of the method.
Although the pinning parameters restrict the system from making certain
transitions, which typically slows down the dynamics, especially at low
temperatures, they are necessary to avoid a certain kind of oscillatory
behavior characterized by an alternation between two distinct spin
configurations as the temperature drops close to zero, see [7, Fig. 1]. Such
oscillatory behavior tends to manifest when the pinning parameters assume
values between $0$ and $\lambda/2$. This phenomenon can be widely observed,
for example, in anti-ferromagnetic Ising models, due to the fact that, in this
case, each spin tends to misalign with its neighbors, in such a way that, if
the pinning parameters are not set appropriately, it may result in an
oscillation where the majority of spins are simultaneously updated to $+1$,
and in the next step, the majority of the spins are updated to $-1$, and so
on. Similarly, for the $\varepsilon$-SCA, if the parameter $\varepsilon$ takes
values in an interval close to $1$, more spin-flips per update are allowed,
and its behavior tends to be similar to the behavior of an SCA without pinning
parameters, which in some cases may also introduce the same oscillatory
behavior. For that reason, more theoretical studies are necessary in order to
derive an analogous assumption for the parameter $\varepsilon$ so that we can
have a theoretical limit for its value below which the convergence is assured.

##  3 Simulations

In this section we exemplify the application of the algorithms presented in
Section 1 in the search for ground states corresponding to three famous
problems:

  1. 1.

_Spin glasses._ Let us consider a spin glass Hamiltonian in a complete graph
with $N$ vertices, without external fields (i.e., $h_{x}=0$), where the values
for the spin-spin couplings $J_{x,y}=J_{y,x}$ are realizations of i.i.d.
random variables. Two cases are analyzed: the Gaussian case, where each
$J_{x,y}$ is a realization of a standard normal random variable, and the
Bernoulli case, where $\mathbb{P}(J_{x,y}=+1)=1-\mathbb{P}(J_{x,y}=-1)=p$.

  2. 2.

_Max-cut problem._ The Hamiltonian is defined in an Erdös-Rényi random graph
with $N$ vertices and edge probability $p$, without external fields, and spin-
spin coupling satisfying $J_{x,y}=-1$ if $\\{x,y\\}$ is an edge of the graph
and $J_{x,y}=0$ otherwise.

  3. 3.

_Traveling salesman problem (TSP)._ The Hamiltonian considered corresponds to
the problem of finding the shortest path that connects $n$ cities, visiting
each city only once, and returning to the initial point, see [16] for more
details.

Despite the fact one of our theoretical results (Theorem 2.3) and some of the
well-established theorems on simulated annealing (e.g., [12]) ensure the
convergence of certain annealing algorithms to the ground states of a
Hamiltonian when the temperature decreases at a logarithmic rate towards zero,
practical simulations involving parallel spin updates have been successfully
applied by decreasing the temperature exponentially fast, e.g., [7, 17, 21].

In the following, we consider the problems described above in the particular
case where the underlying graph contains $N=100$ vertices and compare the
performances of simulated annealing algorithms based on the $\varepsilon$-SCA,
SCA (with pinning parameters $q_{x}=\lambda/2$) and Glauber dynamics. For each
comparison, we ran $1000$ trials for each algorithm, where each trial
consisted of taking an initial spin configuration uniformly at random and
applying $10^{4}$ Markov chain steps with the exponential temperature cooling
schedule given by

$\beta_{t}=\beta_{0}\exp{(\alpha t)},$ (3.1)

where $\beta_{0}=\alpha=10^{-3}$. After that, the histogram of the smallest
energy reached within each trial was generated. In the particular case where
we analyze the TSP, for the sake of simplicity, we considered only $n=10$
cities, that is, $N=n^{2}=100$. The distances between cities $i$ and $j$,
namely $d_{i,j}=d_{j,i}$, were randomly selected integers between $1$ and
$100$, and the parameters $A$ and $B$ (see [16]) were chosen as
$A=\max_{i,j}d_{i,j}$ and $B=1$ .

The histograms respective to each studied case are illustrated in Figs. 1 and
2, and results of the simulations are summarized in Table 1.

(a) Gaussian case, where $\varepsilon=0.9$

(b) Bernoulli case with $p=0.2$, where $\varepsilon=0.35$

(c) Bernoulli with $p=0.5$, where $\varepsilon=0.75$

Fig. 1: Histograms obtained by using the $\varepsilon$-SCA, SCA, and Glauber
dynamics for spin glasses on the complete graph with $N=100$.

(a) Max-cut problem with $p=0.1$, where $\varepsilon=0.6$

(b) Max-cut problem with $p=0.9$, where $\varepsilon=0.1$

(c) TSP, where $\varepsilon=0.3$

Fig. 2: Histograms obtained through the $\varepsilon$-SCA, SCA, and Glauber
dynamics applied to max-cut problems on Erdös-Rényi random graphs and the TSP,
where $N=100$. Table 1: Success rate of obtaining a ground state associated to
different Hamiltonians by using $\varepsilon$-SCA, SCA, and Glauber dynamics

Models | $\min(H)$ |  $\varepsilon$-SCA | SCA |  | Glauber
---
dynamics
| Gaussian
---
spin glass
$-739.5749$ | $89.5\%$ | $31.5\%$ | $5.5\%$
| Bernoulli
---
spin glass
($p=0.2$)
$-600$ | $88.3\%$ | $0\%$ | $5.5\%$
| Bernoulli
---
spin glass
($p=0.5$)
$-746$ | $42.9\%$ | $22.6\%$ | $3.6\%$
| Bernoulli
---
spin glass
($p=0.8$)
$-2944$ | $100\%$ | $100\%$ | $100\%$
| Max-cut
---
problem
($p=0.1$)
$-219$ | $73.3\%$ | $11.6\%$ | $5\%$
| Max-cut
---
problem
($p=0.9$)
$-279$ | $52.2\%$ | $0\%$ | $8.7\%$
TSP | -47453 | $3.8\%$ | $0\%$ | $0.1\%$

##  4 Discussion

Through the inspection of Figs. 1 and 2 and the analysis of Table 1, we
immediately observe that the simulated annealing algorithm based on
$\varepsilon$-SCA is the one which offered the highest success rates in
obtaining ground states in all considered scenarios, being consistent with our
previous findings in [6, 7]. Since the study of Gaussian spin glasses and the
max-cut problem on Erdös-Rényi random graphs was already addressed in [7], let
us solely concentrate on the analysis of Bernoulli spin glasses and the TSP.

In the minimization problem of Bernoulli spin glass Hamiltonians and in the
TSP, we verified a phenomenon that is similar to one previously observed in
the max-cut problem on Erdös-Rényi random graphs, where we noticed that the
SCA struggles in obtaining low energy states when antiferromagnetic spin-spin
interactions (i.e., when $J_{x,y}<0$) are more prevalent than ferromagnetic
interactions (i.e., when $J_{x,y}>0$). In these cases, the pinning parameters,
homogeneously taken as $q_{x}=\lambda/2$, prevent the system from reaching
certain low-energy configurations. Therefore, in order to overcome such a
limitation, finding optimal values for the pinning parameters by refining the
results from [17] or increasing the simulation times by letting the
temperature decrease at a slower rate can serve as alternatives for improving
the accuracy of the SCA. In the cases where such a phenomenon described above
was not present, the SCA reached ground states at a higher rate compared to
the Glauber dynamics.

The $\varepsilon$-SCA was applied to Bernoulli spin glasses with probabilities
$p$ chosen as $p=0.2,0.5$ and $0.8$, considering the values for $\varepsilon$
equal $\varepsilon=0.35,0.75$ and $1$, respectively. For the case $p=0.2$
(resp. $p=0.8$), where antiferromagnetic (resp. ferromagnetic) interactions
are dominant, the $\varepsilon$-SCA obtained a higher success rate in finding
a ground state of $88.3\%$ (resp. $100\%$). On the other hand, for the case
$p=0.5$, where the number of ferromagnetic and antiferromagnetic interactions
are expected to be the same, the performance of the $\varepsilon$-SCA was not
as high as in the other two cases. For such a specific model, a ground state
seems to be a spin configuration whose number of $+1$ and $-1$ spins are
nearly the same, whereas there is a large number of configurations sharing
that same feature whose energy is close to the energy of the ground state. So,
for some reason related to the dispositions of the ferromagnetic and
antiferromagnetic interactions, there is an increased chance for the dynamics
of getting trapped in one of these configurations and not converging to the
ground state. Again, considering slower decreasing temperature schedules may
also improve the accuracy of the algorithm.

Although the success rate obtained for the $\varepsilon$-SCA when applied to
solve the TSP turned out to be smaller compared to the obtained when applied
to the other problems, this value is still greater than those obtained when
SCA and Glauber dynamics were applied instead. The $\varepsilon$-SCA, with
$\varepsilon=0.3$, and the Glauber dynamics obtained the same lowest energy
configuration whose energy was equal $-47453$, corresponding to a trajectory
with length equal $97$, with success rates of $3.8\%$ and $0.1\%$,
respectively. On the other hand, the lowest energy configuration obtained for
the SCA had energy equal to $-47369$, which corresponded to a trajectory with
a length equal to $181$. Since the TSP can be expressed as the minimization
problem of an Ising Hamiltonian with antiferromagnetic interactions and
negative external fields, its energy landscape has been revealed to be much
more complex compared to the other Hamiltonians. Such fact leads the algorithm
to get stuck in configurations that are not ground states if we do not allow
the temperature to drop at a sufficiently low speed. Moreover, the choice of
the parameters $A$ and $B$ from its Hamiltonian may also affect the
performance of the algorithms since a large value for the ratio $A/B$ may
increase the chances of the dynamics stopping in a configuration that
corresponds to a possible trajectory for the salesman but with non-minimal
length, however, no general method for determining their optimal values is
known. As an attempt to mitigate the difficulty of reaching a ground state at
a higher rate without having to increase simulation times, we compared the
performance of the $\varepsilon$-SCA with the so-called Digital Annealer’s
Algorithm (often referred to as DA) [1] applied to the TSP. Such a comparison
was performed under the same conditions as before, where the latter algorithm
obtained the same ground state with a success rate of $8\%$, see Fig. 3. In
that case, the DA obtained not only a ground state at a higher rate, but also
obtained low energy states more frequently. It follows that we found a
particular problem where $\varepsilon$-SCA’s performance can be surpassed by a
single spin-flip algorithm. Even though the DA offers superior performance,
there is still no theoretical explanation for the reason why this algorithm
performs better, specifically in this kind of problem.

Fig. 3: Histograms corresponding to the application of the $\varepsilon$-SCA
(with $\varepsilon=0.3$) and the DA to the TSP.

Similarly as in [7], let us address the topic regarding the optimal value of
$\varepsilon$ for the $\varepsilon$-SCA that should be chosen in order to
maximize the success rate of reaching a ground state. The simulations were
performed by considering the same Bernoulli spin glass Hamiltonians, with the
same cooling schedule and simulation time as those used previously. In the
plots from Fig. 4, each point that corresponds to a certain value of
$\varepsilon$ represents the success rate for the $\varepsilon$-SCA in
obtaining a ground state within $1000$ trials. Let us mention that in the case
$p=0.8$ the algorithm converged to the ground state with a success rate of
$100\%$ for all tested values of $\varepsilon$ in the interval $(0,1]$, for
that reason, we did not illustrate the plot corresponding to this case. On the
other hand, in Figs. 4(a) and 4(b), we observe a growth tendency of such a
rate as $\varepsilon$ increases, but an abrupt decrease appears as
$\varepsilon$ assumes values in a region close to $1$, which corresponds to
the region where oscillation occurs. Note that for the Bernoulli spin glass
with $p=0.2$ such decrease starts occurring for smaller values of
$\varepsilon$ in comparison with the case $p=0.5$ due to the fact that anti-
ferromagnetic interactions dominate for smaller values of $p$ and the system
becomes more susceptible to oscillate at low temperatures provided a larger
number of spin-flips are allowed. Furthermore, as we consider larger values
for $p$ so that ferromagnetic interactions play a more significant role in the
Hamiltonian, the ground state tends to approach the configuration whose spins
are all the same (either all $+1$ or all $-1$), and, in that case, the
algorithm will converge to the ground states independently on the value of
$\varepsilon$ in the interval $(0,1]$.

(a) $N=100$ and $p=0.2$

(b) $N=100$ and $p=0.5$

Fig. 4: The dependence on $\varepsilon$ of the success rate of finding a
ground state for the Bernoulli spin glass Hamiltonian.

##  5 Conclusion

In this work we extended the preliminary results showed in [7] by providing
rigorous and practical results that can serve as a justification for the
employment of simulated annealing algorithms based on parallel spin-flip
updates for obtaining ground states. As a theoretical addition, it was
possible to show that the mixing time for the $\varepsilon$-SCA as well as for
the SCA is at most proportional to $\log|V|$ provided the temperature is
sufficiently high.

The applications of simulated annealing for the $\varepsilon$-SCA, SCA, and
Glauber dynamics were illustrated by considering the minimization problems of
Gaussian and Bernoulli spin glass Hamiltonians, the max-cut problem on Erdös-
Rényi random graphs, and the traveling salesman problem (TSP). The SCA
outperformed the Glauber dynamics in some cases but failed in obtaining ground
states for certain models where antiferromagnetic interactions play a
significant role such as in max-cut problems on Erdös-Rényi random graphs with
values for $p$ in an interval close to 1, Bernoulli spin glasses with $p$ in
an interval close to $0$, and the TSP. However, the $\varepsilon$-SCA
exhibited the largest success rates in all scenarios, being consistent with
our previous findings [6, 7].

Even though the $\varepsilon$-SCA performed significantly better than the
other two algorithms in all cases, its hypothetical advantage of allowing a
large number of spin flips and not suffering from the effect of pinning
parameters did not guarantee a large success rate for the TSP. Due to its
complexity, in order to improve the results for the TSP, larger simulation
times with slowly decreasing temperatures are required. Moreover, for this
specific class of problems, the so-called Digital Annealer’s Algorithm (DA)
showed to be more efficient than $\varepsilon$-SCA. For that reason,
theoretical and practical aspects concerning the properties of the DA are
worth investigating.

The next step in future investigations consists of theoretically justifying
the efficiency and limitations of the use of exponential cooling schedules for
the SCA as well as for the $\varepsilon$-SCA, aiming at providing some
theoretical background that can be decisive in real-world applications. In
order to gain more information that may allow us to classify the algorithms
according to their appropriateness in solving certain classes of problems,
more simulations considering a larger variety of problems including larger
values of $N$ will be necessary. In order to do so, the employment of hardware
accelerators such as GPUs are being considered, and a brief evaluation has
been reported in [6].

## References

  * [1] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, H.G. Katzgraber.  Physics-Inspired Optimization for Quadratic Unconstrained Problems Using a Digital Annealer.  _Frontiers in Physics_ , 7 (2019).
  * [2] F. Barahona, M. Grötschel, M. Jünger, G. Reinelt.  An Application of Combinatorial Optimization to Statistical Physics and Circuit Layout Design.  _Operations Research_ , 36 (1988) 493–513.
  * [3] O. Catoni.  Rough Large Deviation Estimates for Simulated Annealing: Application to Exponential Schedules.  _The Annals of Probability_ , 20 (1992): 1109–1146.
  * [4] V. Černý.  Thermodynamical approach to the traveling salesman problem: an efficient simulation algorithm.  _J. Optimiz. Theory App._ , 45 (1985): 41–51.
  * [5] J. Ding and Y. Peres.  Mixing time for the Ising model: A uniform lower bound for all graphs.  _Annales de l’Institut Henri Poincaré, Probabilités et Statistiques_ , 47 (2011) 1020–1028.
  * [6] B.H. Fukushima-Kimura, S. Handa, K. Kamakura, Y. Kamijima, K. Kawamura, A. Sakai.  Mixing time and simulated annealing for the stochastic cellular automata.  arXiv preprint arXiv:2007.11287.
  * [7] B.H. Fukushima-Kimura, Y. Kamijima, K. Kawamura and A. Sakai.  Stochastic optimization via parallel dynamics: rigorous results and simulations.  _Proceedings of the ISCIE International Symposium on Stochastic Systems Theory and its Applications_ , 2022, 65-71.
  * [8] M.R. Garey and D.S. Johnson.  Computers and Intractability: A Guide to the Theory of NP-Completeness.  W.H. Freeman and Company, San Francisco (1979).
  * [9] R.J. Glauber.  Time-dependent statistics of the Ising Model.  _J. Math. Phys._ , 4 (1963): 294–307.
  * [10] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, K. Tatsumura.  High-performance combinatorial optimization based on classical mechanics.  _Science Advances_ , 7 (2021) eabe7953.
  * [11] H. Goto, K. Tatsumura, A.R. Dixon.  Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems.  _Science Advances_ , 5 (2019) eaav2372.
  * [12] B. Hajek.  Cooling schedules for optimal annealing.  _Math. Oper. Res._ , 13 (1988): 191–376.
  * [13] S. Handa, K. Kamakura, Y. Kamijima and A. Sakai.  Finding optimal solutions by stochastic cellular automata.  Preprint. arXiv:1906.06645.
  * [14] T.P. Hayes and A. Sinclair  A general lower bound for mixing of single-site dynamics on graphs.  _The Annals of Applied Probability_ , 17 (2007) 931–952.
  * [15] D.A. Levin and Y. Peres.  Markov Chains and Mixing Times, second edition.  AMS, Providence, Rhode Island (2017).
  * [16] A. Lucas.  Ising formulations of many NP problems.  _Front. Phys._ , 12 (2014): https://doi.org/10.3389/fphy.2014.00005.
  * [17] T. Okuyama, T. Sonobe, K. Kawarabayashi, and M. Yamaoka.  Binary optimization by momentum annealing.  _Phys. Rev. E_ , 100 (2019), 012111 .
  * [18] P. Dai Pra, B. Scoppola and E. Scoppola.  Sampling from a Gibbs measure with pair interaction by means of PCA.  _J. Stat. Phys._ , 149 (2012): 722–737.
  * [19] H. Sakaguchi, K. Ogata, T. Isomura, S. Utsunomiya, Y. Yamamoto, K. Aihara.  Boltzmann sampling by degenerate optical parametric oscillator network for structure-based virtual screening.  _Entropy_ , 18 (2016).
  * [20] B. Scoppola and A. Troiani.  Gaussian mean field lattice gas.  _J. Stat. Phys._ , 170 (2018): 1161–1176.
  * [21] K. Yamamoto, K. Kawamura, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, H. Teramoto, A. Sakai, S. Takamaeda-Yamazaki, M. Motomura.  STATICA: A 512-Spin 0.25M-Weight Annealing Processor With an All-Spin-Updates-at-Once Architecture for Combinatorial Optimization With Complete Spin–Spin Interactions.  _IEEE Journal of Solid-State Circuits_ , 56 (2021): 165–178.