Simulated bifurcation assisted by thermal fluctuation

Unconventional computing with special-purpose hardware devices for solving combinatorial optimization problems has attracted growing interest due to practical importance. Many combinatorial optimization problems can be mapped onto finding ground states of Ising spin models Barahona1982 ; Lucas2014 , which is referred to as the Ising problem. Special-purpose hardware devices for the Ising problem are called Ising machines. Ising machines utilizing natural phenomena have been developed, such as quantum annealers Kadowaki1998 ; Farhi2000 ; Farhi2001 ; Das2008 ; Albash2018 ; Hauke2020 with superconducting circuits Johnson2011 , coherent Ising machines with pulse lasers Wang2013 ; Marandi2014 ; Inagaki2016 ; Yamamoto2017 ; Honjo2021 , oscillator-based Ising machines Wang2017 ; Wang2019a ; Chou2019 ; Mallick2020 ; Vaidya2022 , and other Ising machines with various systems, such as stochastic nanomagnets Sutton2017 , gain-dissipative systems Kalinin2018 , spatial light modulators Pierangeli2019 , memristor Hopfield neural networks Cai2020 , and spin torque nano-oscillators Houshang2020 ; Albertsson2021 .

Ising machines have also been implemented with special-purpose digital processors Yamaoka2016 ; Tsukamoto2017 ; Aramon2019 ; Okuyama2019 ; Yamamoto2021 ; Patel2020 ; Leleu2021a using simulated annealing (SA) Kirkpatrick1983 and other algorithms. Among such algorithms, simulated bifurcation (SB) is a recently proposed heuristic algorithm Goto2019 . SB originates from numerical simulations of Hamiltonian dynamics with bifurcation that is a classical counterpart of quantum adiabatic bifurcation in nonlinear oscillators Goto2016 ; Goto2019a , which itself has been studied actively Nigg2017 ; Puri2017 ; Zhao2018 ; Goto2018a ; Kewming2020 ; Onodera2020 ; Goto2020b ; Kanao2021 ; Goto2021b . Simulation-based approaches such as SB allow one to deal with dense spin-spin interactions with high precision, which might be challenging for physical implementations. Also, the classical dynamics can be simulated efficiently, unlike quantum dynamics. SB can be accelerated by parallel processing with, e.g., field-programmable gate arrays (FPGAs) Goto2019 ; Tatsumura2019 ; Zou2020 ; Tatsumura2021 , because of its capability of simultaneous updating of variables. Recently proposed variants of SB have achieved faster and more accurate optimization Goto2021 than original SB.

To further improve the performance of SB, here we introduce thermal fluctuation to SB. SA can yield high-accuracy solutions by modeling thermal fluctuation Kirkpatrick1983 , while quantum annealing utilizes quantum fluctuation Kadowaki1998 . These fluctuations can assist escape from local minima of the Ising problem and lead to higher solution accuracy.

Our method is based on the Nosé-Hoover method Nose1984 ; Hoover1985 , which enables simulations of Hamiltonian dynamics at finite temperatures Leimkuhler2004 . Unlike SA, the Nosé-Hoover method does not use random numbers, namely, is deterministic, and thus the simplicity of SB is preserved. We find that a simplified dynamics with only a heating process can improve the performance of SB, where an ancillary dynamical variable in the Nosé-Hoover method is replaced by a constant. We numerically demonstrate this improvement by solving instances of the Ising problem with up to 2000 spin variables and all-to-all connectivity, which corresponds to the Sherrington-Kirkpatrick (SK) model introduced in studies of spin glasses Sherrington1975 ; Parisi1979 ; Parisi1993 . The SK model has been widely used to measure the performance of Ising machines Inagaki2016 ; Aramon2019 ; Honjo2021 ; Goto2019 ; Okuyama2019 ; Goto2021 ; Yamamoto2021 ; Leleu2021a ; Oshiyama2022 . Proposed heated SB is also suitable for massively parallel implementations with, e.g., FPGAs.

Here we apply the Nosé-Hoover method Nose1984 ; Hoover1985 ; Leimkuhler2004 with a finite temperature $T$ to Eqs. (2) and (3), obtaining

where $\xi$ is an ancillary variable playing a role of thermal fluctuation, and $M$ is a parameter (mass). The variable $\xi$ controls an instantaneous temperature defined by

to be a given $T$ as follows. When $T_{\rm inst}$ is smaller than $T$, $\dot{\xi}$ is negative according to Eq. (8), which makes $\xi$ negative. Then ${|y_{i}|}$ increase owing to the last term in Eq. (7), and thus $T_{\rm inst}$ increases and approaches $T$, which can be regarded as heating. To the contrary, when ${T_{\rm inst}>T}$, cooling occurs.

We found that SB gives better solutions when $\xi$ is kept negative by negative initial $\xi$ and large $M$ (leading to ${\dot{\xi}\simeq 0}$). This observation suggests that the heating can improve SB but the cooling is unnecessary. This may be because increased ${\left|y_{i}\right|}$ by the heating can lead to escape from local minima of the Ising problem. Furthermore, small ${|\dot{\xi}|}$ due to large $M$ above implies that constant $\xi$ can play a similar role, and then we found that $\xi$ replaced by a negative constant $-\gamma$ ${(\gamma>0)}$ can rather yield higher performance. The constant $\gamma$ is regarded as a rate of the heating.

Thus, in this paper, we propose SB with a heating term, which we call heated bSB (HbSB) and dSB (HdSB), as follows,

We numerically solve Eqs. (10) and (11) by discretizing the time by ${t_{k+1}=t_{k}+\Delta t}$ with a time interval $\Delta t$, and by calculating $x_{i}\left(t_{k+1}\right)$ and $y_{i}\left(t_{k+1}\right)$ from $x_{i}\left(t_{k}\right)$ and $y_{i}\left(t_{k}\right)$ in each time step. Here note that the symplectic Euler method is not applicable for Eqs. (10) and (11), because these equations are no longer Hamiltonian equations owing to the term $\gamma y_{i}$ Leimkuhler2004 . We empirically found that solution accuracy can be improved by the same update as previous bSB and dSB Goto2021 followed by an update corresponding to the term $\gamma y_{i}$. This ordering results in nonzero momenta by the heating and can prevent from getting stuck at the walls. See Methods for a detailed algorithm.

In the following, we compare heated SB with previous SB by solving instances of the Ising problem with all-to-all connectivity, where $J_{ij}$ are randomly chosen from $\pm 1$ with equal probabilities (corresponding to the SK model). The control parameter $a(t)$ is linearly increased from 0 to $a_{0}$. The constant parameters are set as ${a_{0}=1}$ and

where $c_{1}$ is a parameter tuned around 0.5, which is based on random matrix theory Goto2019 ; Goto2021 . $x_{i}$ and $y_{i}$ are initialized by uniform random numbers in the interval ${(-1,1)}$.

To confirm the effect of the heating, we calculate the instantaneous temperature $T_{\rm inst}$ in Eq. (9) and the cut value $C$ in Eq. (14) at every time step. Figure 1a shows typical examples of time evolutions of $T_{\rm inst}$. Here parameters are the ones optimized in advance, which are explained later. For bSB, $T_{\rm inst}$ rapidly decreases owing to collisions with the perfectly inelastic walls, while $T_{\rm inst}$ for HbSB is kept much higher owing to the heating, as expected. For dSB, $T_{\rm inst}$ is also higher than bSB, because the discretized forces $f_{i}$ in Eq. (5) increase energy, violating conservation of energy Goto2021 . In comparison between HbSB and dSB, $T_{\rm inst}$ for HbSB is higher than dSB around the end of the evolution. HdSB shows similar $T_{\rm inst}$ to dSB, because the optimal rate of the heating $\gamma$ for HdSB is small.

Figure 1b shows last parts of typical time evolutions of $C$. For bSB, $C$ is almost constant for the last 1000 time steps, while $C$ for HbSB continues to fluctuate until nearly the end of the evolution owing to the heating. The fluctuation for HbSB looks similar to the fluctuations for dSB and HdSB.

Next, we solve K₂₀₀₀ in 10⁴ trials with random initial $x_{i}$ and $y_{i}$. In one trial, $C$ is evaluated every 100 time steps, and the best value is output. Figure 2 shows examples of cumulative distributions of $C$, where the number of trials giving cut values lower than $C$ are normalized by the total number of trials, 10⁴. For bSB, the majority of trials results in one of a few values of $C$ between 33200 and 33300. On the other hand, dSB, HbSB, and HdSB yield broader distributions, owing to fluctuations (higher instantaneous temperatures) in their dynamics. It is notable that the heating makes the distribution of bSB similar to dSB (and HdSB). Inset in Fig. 2 shows a magnification around the best known cut value, ${C_{\rm best}=33337}$ Goto2021 . The distributions for HbSB and HdSB are shifted to larger $C$ than dSB, indicating improvement by the heating.

We then evaluate average and maximum $C$ in the $10^{4}$ trials and probability $P$ for obtaining $C_{\rm best}$. Here, $P$ is estimated by dividing the number of trials obtaining $C_{\rm best}$ by the total number of trials. Also, using $P$ and the number of time steps, $N_{\rm s}$, we calculate the number of time steps required to find $C_{\rm best}$ with a probability of 99%, which we call step-to-solution $S$, given by

Step-to-solution is a useful measure of performance of an algorithm Leleu2021a . (Time-to-solution is often used to measure performance of Ising machines Aramon2019 ; Goto2021 ; Leleu2021a ; Patel2020 , but it depends on not only algorithms but also implementations to hardware devices. Time-to-solution equals step-to-solution multiplied by a computation time for one time step.) Here the parameters $\Delta t$, $c_{1}$, and $\gamma$ are set such that $S$ are minimized. For bSB, instead, average $C$ is maximized, because we found $P=0$ for bSB and could not estimate $S$.

Figure 3a shows average and maximum $C$ as functions of $N_{\rm s}$. For large $N_{\rm s}$, average $C$ for dSB, HbSB, and HdSB are larger than that for bSB. $C_{\rm best}$ is reached by three SBs other than bSB. Figure 3b shows that, for large $N_{\rm s}$, $P$ are the highest for HbSB, followed by HdSB, dSB, and bSB in this ordering.

Figure 3c shows $S$ as functions of $N_{\rm s}$. Each SB has a minimum of $S$ at certain $N_{\rm s}$, and in the following we compare $S$ minimized with respect to $N_{\rm s}$. The black cross represents a value obtained from previously reported data for dSB Goto2021 . The present value of $S$ for dSB is smaller than the previous value, because in this study the parameters $\Delta t$ and $c_{1}$ are optimized for K₂₀₀₀ while not in the previous study. In Fig. 3c, at optimal $N_{\rm s}$, $S$ for HbSB is the smallest. Compared with dSB, $S$ are reduced by 32.7% for HbSB and 19.7% for HdSB. These results demonstrate that the heating improves the performance. Although dSB performs better than bSB, bSB is much more improved by the heating than dSB, and resulting HbSB shows higher performance than HdSB.

Figure 4 shows the medians of $S$ for the 100 instances Goto2021 ; Leleu2021a as functions of $N_{\rm s}$. HbSB results in the smallest $S$ among four SBs ($S$ by bSB is much larger than those by the others). In comparison with dSB, HbSB reduces $S$ by 9.55%. This result demonstrate that HbSB can improve the performance for not only K₂₀₀₀ but also the other instances.

We have demonstrated that SB can be improved by introducing fluctuation with a heating term, which has been obtained by replacing an ancillary dynamical variable in the Nosé-Hoover method by a constant rate of heating. This heating can be effective for both bSB and dSB, but the improvements are much larger for bSB. We have solved all-to-all connected 2000-spin and 700-spin instances of the Ising problem (the SK model) and have found that HbSB gives better step-to-solution than bSB, dSB, and HdSB. Since proposed heated SB shares the simple dynamics with previous SB, we expect that heated SB will be accelerated by massively parallel processing implemented by, e.g., FPGAs. This study also implies that further improvements of SB will be possible by simple physics-inspired modifications like the heating term introduced here.