Research ArticlePHYSICS

High-performance combinatorial optimization based on classical mechanics

See allHide authors and affiliations

Science Advances  03 Feb 2021:
Vol. 7, no. 6, eabe7953
DOI: 10.1126/sciadv.abe7953


Quickly obtaining optimal solutions of combinatorial optimization problems has tremendous value but is extremely difficult. Thus, various kinds of machines specially designed for combinatorial optimization have recently been proposed and developed. Toward the realization of higher-performance machines, here, we propose an algorithm based on classical mechanics, which is obtained by modifying a previously proposed algorithm called simulated bifurcation. Our proposed algorithm allows us to achieve not only high speed by parallel computing but also high solution accuracy for problems with up to one million binary variables. Benchmarking shows that our machine based on the algorithm achieves high performance compared to recently developed machines, including a quantum annealer using a superconducting circuit, a coherent Ising machine using a laser, and digital processors based on various algorithms. Thus, high-performance combinatorial optimization is realized by massively parallel implementations of the proposed algorithm based on classical mechanics.


Combinatorial optimization problems appear in various social and industrial situations, so quickly solving such problems makes the society and industry more efficient. However, these problems are notoriously hard due to combinatorial explosion, an exponential increase in the number of candidate solutions depending on the problem size (1). Thus, novel computational approaches to combinatorial optimization have been expected. A well-known example is a quantum annealer (QA), which is based on quantum annealing (24) and its superconducting circuit implementation (5, 6). The QA is an Ising machine designed to find ground states of Ising spin models (7). Such Ising machines are believed to be broadly useful, because the Ising problem belongs to the nondeterministic polynomial time (NP)–complete class (7), and consequently, many combinatorial optimization problems can be reduced to the Ising problem (8). Various Ising machines other than QA have been developed, including coherent Ising machines (CIMs) implemented with pulse lasers (914) and other kinds of optical Ising machine (1518), as well as digital processors based on simulated annealing (SA) (1923), which is a standard heuristic algorithm for combinatorial optimization (1, 24, 25), or other recently proposed algorithms (2629).

A heuristic algorithm called simulated bifurcation (SB) has recently been proposed for accelerating combinatorial optimization (30). SB is a purely quantum-inspired algorithm, that is, it was derived from a classical-mechanical model corresponding to a quantum computer called a quantum bifurcation machine (QbM) (3033), which is based on quantum adiabatic optimization using nonlinear oscillators exhibiting quantum-mechanical bifurcation phenomena. Consequently, SB is based on numerical simulation of adiabatic evolutions in classical nonlinear Hamiltonian systems exhibiting bifurcations (34). Different dynamical approaches have also recently been proposed (3541). SB is not based on the gradient method, unlike other dynamical approaches such as the Hopfield-Tank model (42), simulated CIM (SimCIM) (35, 39), and their variants, but based on adiabatic evolutions of energy conservative systems like purely adiabatic QA and QbM. [Such interesting contrast between QbM and CIM has been summarized in a review paper on this topic (33).]

Unlike SA in general cases, SB allows simultaneous updating of variables and therefore can easily accelerate combinatorial optimization through massively parallel processing using modern many-core processors such as field-programmable gate arrays (FPGAs) (30, 43, 44) and graphic processing units (GPUs) (30). An SB-based machine (SBM) implemented with a single FPGA, where about 8000 operations are performed in parallel, was able to find good approximate solutions of a 2000-spin Ising problem in 0.5 ms, about 10 times faster than a CIM (30). This result suggests that parallelizability is a key property of optimization algorithms for their acceleration by fully exploiting modern high-performance computing systems. In this direction of research, other parallelizable algorithms have also recently been proposed by mapping a given problem to a bipartite one and applying parallel SA updating to each group of spins (2628). These new algorithms essentially rely on the same mechanism, although they are given different names: momentum annealing (MA) (26), stochastic cellular automata annealing (SCA) (27), and restricted Boltzmann machine (RBM)’s parallel stochastic sampling (28).

The previous results on SB demonstrate that SB is useful for quickly finding good approximate solutions. However, it remains unclear whether SB can find optimal solutions of large-scale problems. For enhancing the power of SB in terms of solution accuracy, in this work, we introduce two SB variants, named ballistic SB (bSB) and discrete SB (dSB), in addition to the original adiabatic SB (aSB). We solve various problems to compare the performance of bSB and dSB with that of aSB and other recently developed machines, including a QA, a CIM, and digital processors based on various algorithms. This benchmarking shows that bSB and dSB provide faster and more accurate optimizations than does aSB and that our new SBMs achieve high performance compared to the other machines. dSB can find optimal or near-optimal solutions of problems with up to one million spins, which aSB and bSB cannot achieve. Thus, high-performance machines for combinatorial optimization are realized by massively parallel implementations of the proposed algorithms based on classical mechanics.


bSB and dSB algorithms

The Ising problem is to find a spin configuration minimizing the Ising energy, defined asEIsing=12i=1Nj=1NJi,jsisji=1Nhisi(1)where si denotes the ith spin taking 1 or −1, N is the number of spins, Ji,j is the coupling coefficient between the ith and jth spins (Ji,j = Jj,i and Ji,i = 0), and hi is the local field on the ith spin. Since introducing an ancillary spin reduces the Ising problem to the one without local fields (see section S1), here, we focus on the Ising problem with no local fields (hi = 0).

To solve the Ising problem, QbM uses quantum-mechanical nonlinear oscillators called Kerr-nonlinear parametric oscillators (KPOs), each of which can generate a Schrödinger cat state, i.e., a quantum superposition of two oscillating states, via a quantum-mechanical bifurcation (31). Such a KPO has recently been realized experimentally using superconducting circuits (45, 46). However, the realization of a large-scale QbM will require a long time. On the other hand, it was found that a classical-mechanical model corresponding to QbM, referred to as a classical bifurcation machine (CbM), also works as an Ising machine (31, 33). In this case, we can efficiently simulate such a classical machine using present digital computers, instead of building real machines. [This is not the case for QbM, because QbM can also be used as a universal quantum computer (47, 48), which is classically unsimulatable.] This simulation approach paves the way for large-scale Ising machines with all-to-all connectivity.

By modifying the equations of motion for CbM such that computational costs are reduced and the symplectic Euler method (49) is applicable, we obtain the following Hamiltonian equations of motion for aSB (30)x·i=HaSByi=a0yi(2)y·i=HaSBxi=[xi2+a0a(t)]xi+c0j=1NJi,jxj(3)HaSB=a02i=1Nyi2+VaSB(4)VaSB=i=1N(xi44+a0a(t)2xi2)c02i=1Nj=1NJi,jxixj(5)where xi and yi are, respectively, the position and momentum of a particle corresponding to the ith spin, dots denote time derivatives, a(t) is a control parameter increased from zero, a0 and c0 are positive constants, and VaSB is the potential energy in aSB.

To qualitatively explain the operation principle of aSB, we show an example of the dynamics in aSB in Fig. 1 (A and B), where the ferromagnetic two-spin Ising problem (J1,2 = J2,1 = 1) with solutions s1 = s2 = ± 1 is solved as the simplest problem. All positions and momenta are randomly set around zero at the initial time. The initial potential has a single local minimum at the origin (top and middle of Fig. 1B), so the particles circulate around the origin (Fig. 1A and middle of Fig. 1B). When a(t) becomes sufficiently large, bifurcations occur, that is, multiple local minima of the potential appear. Then, the particles adiabatically follow one of the minima. Consequently, each xi bifurcates to a positive or negative value (Fig. 1A and bottom of Fig. 1B). Since two local minima corresponding to the two solutions have lower energies and appear earlier than the other two local minima, the particles successfully find one of the solutions. Last, the sign of xi, sgn(xi), gives the ith spin, si, for the solution of the Ising problem. It has been empirically found that aSB works well for much larger-scale and more complex problems (30).

Fig. 1 Dynamics in three SB algorithms.

Two-spin ferromagnetic Ising problem (J1,2 = J2,1 = 1) is solved by aSB, bSB, and dSB with a(t) = 0.01t, a0 = 1, and c0 = 0.2. (A) Time evolutions of x1 and x2 in aSB. (B) Potential energies VaSB in aSB at the initial (top), intermediate (middle), and final (bottom) times. The vertical dotted line in (A) indicates the intermediate time. Curves in white are trajectories between the initial and intermediate times (middle) and between the intermediate and final times (bottom). Local minima of VaSB are shown by + in red. (C to F) Corresponding results for bSB (C and D) and dSB (E and F). (G to I) Time dependences of potential energy and total energy in aSB (G), bSB (H), and dSB (I). In each panel, the inset magnifies the part in the red dashed square. The dashed arrow in the inset in (I) indicates tunneling-like behavior in dSB. The red bold curve in the middle of (F) shows the trajectory corresponding to this inset.

This aSB relies on the fact that the second term in VaSB is approximately proportional to the Ising energy at the final time (30). In this approximation, analog errors arise from the use of continuous variables (positions) instead of discrete variables (spins). These analog errors in aSB may degrade solution accuracy and result in approximate solutions. Such analog errors in different dynamical approaches have also been discussed (38, 40).

To suppress analog errors, we introduce perfectly inelastic walls at xi = ± 1. That is, at each time, we replace xi with its sign, sgn(xi) = ± 1, and set yi = 0 if ∣xi∣ > 1. These walls force positions to be exactly equal to 1 or −1 when a(t) becomes sufficiently large. Moreover, we drop the fourth-order terms in VaSB, because the inelastic walls can play a role similar to the nonlinear potential walls. We thus obtain the following equationsx·i=HbSByi=a0yi(6)y·i=HbSBxi=[a0a(t)]xi+c0j=1NJi,jxj(7)HbSB=a02i=1Nyi2+VbSB(8)VbSB=a0a(t)2i=1Nxi2c02i=1Nj=1NJi,jxixj when xi1 for all xiotherwise VbSB=(9)This artificial dynamical system is the basis of the bSB. In bSB, we solve Eqs. 6 and 7 by the symplectic Euler method, where time is discretized with a time step Δt, together with the updating for the inelastic walls (see Methods for a detailed algorithm). (If we solve these equations by the standard Euler method, instead of the symplectic Euler method, then solution accuracy becomes lower. See section S2.)

Similar modification to the above walls has been proposed for SimCIM (39), but this algorithm is based on the gradient method, like the Hopfield-Tank model (42), and also uses stochastic processes. In contrast, bSB is based on a classical-mechanical system conserving energy except for the inelastic walls and adiabatic changes of energy and also uses deterministic processes except for initial value setting. As a result, the performance of bSB is quite different from that of SimCIM (see section S3 for the comparison between bSB and SimCIM).

In bSB, it is sufficient to increase a(t) to a0. Then, the final potential has only the second term related to the Ising energy. Consequently, the following condition is satisfied for all i at the final time (see section S4)ΔEi=2sij=1NJi,jsj0(10)where si = sgn (xi) is the sign of xi and ΔEi represents the change in the Ising energy for a flip of si. Note that Eq. 10 is a sufficient condition to show that the spin configuration is a local minimum of the Ising problem. Hence, solutions obtained by bSB are at least local minima in the Ising problem. In contrast, this is not necessarily guaranteed in aSB because of its nonlinear potential terms. (This means that solutions obtained by aSB can sometimes be improved by a naïve local search based on sequential spin flips.) This is another reason why bSB should achieve higher accuracy than aSB. Throughout this work, we linearly increase a(t) from 0 to a0 and set a0 to 1.

Here, we show an example of the bSB dynamics using the same two-spin problem as above. The initial potential has a single local minimum at the origin (top of Fig. 1D) and particles circulate around the origin (Fig. 1C and middle of Fig. 1D), as in aSB. In bSB, however, stable points suddenly jump from the origin to the walls at xi = ± 1, which prevents adiabatic evolution. Instead, particles move toward walls in a ballistic manner (Fig. 1C and bottom of Fig. 1D). This ballistic (nonadiabatic) behavior in bSB leads to fast convergence to a local minimum of VbSB and, consequently, to fast optimization.

For further improvement, we introduce another variant of SB by discretizing xj to sgn(xj) in the second term in Eq. 7x·i=HdSByi=a0yi(11)y·i=HdSBxi=[a0a(t)]xi+c0j=1NJi,jsgn(xj)(12)HdSB=a02i=1Nyi2+VdSB(13)VdSB=a0a(t)2i=1Nxi2c0i=1Nj=1NJi,jxisgn(xj) when xi1 for all xiotherwise VdSB=(14)The discrete version of bSB, namely, dSB (see Methods for a detailed algorithm), is expected to achieve higher accuracy than that of bSB, because analog errors are further suppressed by discretization. We found that this discretization is effective for bSB but not for other algorithms such as SimCIM and aSB (see sections S3 and S5).

Note that the singularity on the boundaries between positive and negative regions has been intentionally neglected. This leads to a violation of conservation of energy across boundaries and, hence, to escape from local minima over potential barriers, as shown below. In this sense, dSB goes beyond naïve algorithms based on classical-mechanical systems conserving energy (except adiabatic change), such as aSB and bSB. We also increase a(t) to a0 for the convergence to a local minimum of the Ising problem at the final time, as in bSB (see section S4).

Figure 1 (E and F) shows an example of the dSB dynamics using the same two-spin problem as above. Unlike aSB and bSB, the particles go back and forth between two local minima through the potential barriers (Fig. 1E and middle of Fig. 1F). This is similar to quantum tunneling, as depicted by the inset in Fig. 1I. This tunneling-like behavior is possible due to the above-mentioned neglect of the singularity on the boundaries; otherwise, the potential walls on the boundaries prevent this tunneling. In contrast, conservation of energy prevents such tunneling in aSB and bSB, as suggested by the insets in Fig. 1 (G and H). Thus, it is expected that this tunneling-like behavior will help dSB to escape local minima of the potential, and hence, dSB will outperform aSB and bSB in terms of solution accuracy.

Note that both bSB and dSB maintain the advantage of aSB over SA, namely, high parallelizability. Therefore, they are expected to realize both high speed and high accuracy simultaneously.

Performance for a 2000-spin Ising problem with all-to-all connectivity

To compare the performance of bSB and dSB with that of aSB, we solved a 2000-spin Ising problem with all-to-all connectivity. This problem was named K2000 and previously solved by aSB (30), a CIM (11), and a recently developed digital chip called STATICA (27), which is based on the above-mentioned SCA. This problem can be regarded as a 2000-node MAX-CUT problem (11, 30); so, here, we evaluate performance using cut values (see section S6 for the definition of the cut value and the relation between MAX-CUT and the Ising problem). The best cut value for K2000 is estimated to be 33,337 (see section S7).

The lines and symbols in Fig. 2A show average and maximum cut values, respectively, for 1000 trials as functions of the number of time steps, Nstep. (Throughout the paper, Nstep denotes the total number of time steps for each trial and the values of cost functions (cut values or Ising energies) as functions of Nstep are final values, not intermediate values, in each trail.) The results clearly show that both bSB and dSB outperform aSB in terms of both speed and accuracy. In addition, only dSB obtained the best value. On the other hand, best values obtained by bSB and aSB become lower for larger Nstep. This result suggests that the best values may be obtained accidentally by nonadiabatic processes in bSB and aSB. For large Nstep, dynamics becomes more adiabatic and the chance to obtain better solutions by nonadiabatic processes may be lost.

Fig. 2 Comparison of performance for a 2000-spin Ising problem with all-to-all connectivity (K2000).

(A) Average (lines) and maximum (symbols) cut values for 1000 trials. Nstep is the number of time steps. The time step Δt is set to 0.5 (aSB) or 1 (bSB and dSB). (See Methods for other settings.) (B) Computation times for our bSB and dSB machines (bSBM and dSBM) implemented with single FPGAs and three previously developed machines: CIM (11), aSB machine (aSBM) (30), and STATICA (27). (The results by STATICA are predicted ones by a simulator. The other results are measured ones using real machines.) The cut values are final ones in each trial, where the computation times of bSBM and dSBM are controlled by Nstep. Lines show average cut values for bSBM (blue) and dSBM (red). Crosses show average cut values for the other three machines. Bars show maximum and minimum cut values. The average, maximum, and minimum values are evaluated by 100 trials. The time step Δt is set to 1.25 for both the bSBM and dSBM. (See Methods for other settings.)

We implemented 2048-spin-size bSB and dSB machines (bSBM and dSBM) using single FPGAs (see section S8 for details) and solved K2000 by them. Figure 2B shows the comparison between our machines and the above three other machines (11, 27, 30), where the lines and the crosses show the average values of our machines and the others, respectively, for 100 trials, and the bars show the maximum and minimum values among the 100 trials. (The bars for our machines are shown at only typical computation times.) Only our dSBM obtained the best value in a short time (2 ms), thereby simultaneously realizing both high speed and high accuracy. Also, our bSBM is remarkably fast, about three times faster than STATICA (27), the previously fastest machine for K2000. Note that the results by STATICA for K2000 are predicted values by a simulator, and a real STATICA chip is still 512-spin size (27). On the other hand, in this work, we have implemented faster real machines.

Benchmarking using time-to-solution and time-to-target

To evaluate the computation speed more quantitatively, here, we introduce two metrics: time-to-solution (TTS) and time-to-target (TTT). TTS is a standard metric for evaluating Ising machine speeds (14, 23, 28, 29), defined as the computation time for finding an optimal or best known value with 99% probability. TTT uses a target value, instead of an optimal value, as a good approximate solution. In this work, we define the target as 99% of the optimal or best known value. TTS and TTT are formulated as Tcom log (1 − 0.99)/ log (1 − PS) (14, 23), where Tcom is the computation time per trial and PS is the success probability for finding the optimal (TTS) or target (TTT) value. PS is estimated from experimental results with many trials. When PS > 0.99, TTS and TTT are defined as Tcom.

In the following, we compare TTS and TTT of our 2048-spin-size bSBM and dSBM with those of other recently developed machines shown in Fig. 3A. Since the bSBM can quickly find good approximate solutions and the dSBM can find optimal solutions of large-scale problems, we use the bSBM and dSBM for evaluations of TTT and TTS, respectively. TTS and TTT of other machines are cited or estimated from the data in the literature (see section S9 for details), because we could not use such machines for the present work. This limits the range of instances that can be used for this benchmarking. Also note that some machines are not the latest ones, as mentioned below.

Fig. 3 TTS and TTT benchmarking.

(A) Machines compared in this benchmarking. Superscripts represent reference numbers for data sources. (B) Comparison of TTTs for two instances of 2000-spin size, K2000 and G22. (C) Comparison of TTSs. The values of the first six problems are medians for 100 (10 for CIM and 20 for QA) random instances. The bar graphs compare all TTTs and TTSs in log scale. For this benchmarking, the time step Δt for bSBM and dSBM is set to the best value for each problem among five values (0.25, 0.5, 0.75, 1, and 1.25). [The same value of Δt is used for 100 random instances of each of the first to sixth problems in (C).] The number of time steps, Nstep, is also optimized for TTT or TTS separately. See Methods and table S1 for the values of Δt and Nstep and other detailed information.

Figure 3B shows the results of TTT. For K2000, the TTT of our bSBM (0.26 ms) is much shorter than those of STATICA (27) (1.50 ms, a predicted value by a simulator) and the CIM (11) (1.1 s). As a 2000-spin-size instance with sparse connectivity, we also solved G22, which is one of the well-known MAX-CUT benchmark instances called G-set and was solved by the CIM (11). For G22, the TTT of our bSBM is two orders of magnitude shorter than that of the CIM. These results demonstrate that our bSBM can find good approximate solutions faster than other recently developed machines of the same spin size. (TTTs of our machines for other G-set instances are provided in table S2.)

Next, we show the results of TTS in Fig. 3C. We start with the same two instances, namely, K2000 and G22. The TTS of our dSBM for them are 1.3 and 2.7 s, respectively. While TTS for K2000 has not been reported so far, TTS for G22 was evaluated with a SimCIM implemented on a FPGA (29), which is based on a recently proposed algorithm called chaotic amplitude control (CAC) (29, 40). The TTS of the SimCIM is estimated at 12 s (see section S9). Thus, our dSBM has achieved a shorter TTS for G22 than that of the state-of-the-art machine. (TTSs of our machines for other G-set instances and the comparison between them and those of the SimCIM are provided in table S2 and fig. S6, respectively.)

We also solved other various instances of the Ising problem (MAX-CUT) by dSBM to compare it with other machines shown in Fig. 3A. For small-scale problems, we can simultaneously perform multiple trials using the 2048-spin-size machine by a block-diagonal structure of the J matrix, as done using a CIM (14). This so-called batch processing improves the success probability PS by selecting the best result among multiple trials, while Tcom is defined as the computation time per batch. In the limit as the number of trials per batch Nbatch goes to infinity, PS may exceed 0.99, and then TTS and TTT become Tcom from the above definitions. In this sense, the TTS and TTT are well defined even for batch processing.

As shown in Fig. 3C, for two small-scale problems with 60 spins (all-to-all connectivity) and 200 spins (sparse connectivity), our dSBM achieved much shorter TTSs than those of a QA and a CIM (14). (This QA is not the latest version.) For 100-spin and 700-spin problems with all-to-all connectivity, the TTSs of the dSBM are much shorter than those of the SimCIM with CAC (29). These short TTSs of dSBM compared to the SimCIM come not from computation speed or implementations but the algorithmic advantage of dSB over CAC. That is, dSB needs fewer matrix-vector multiplications (MVMs), which is the most computation-intensive part in both algorithms, to obtain solutions than does the SimCIM with CAC. [The numbers of MVMs to solutions of the 100-spin and 700-spin problems are 9.4 × 10 and 8.1 × 104, respectively, for dSB but 5.6 × 103 and 7.8 × 105, respectively, for CAC (29).] For two 1024-spin problems (all-to-all connectivity) with different ranges of Ji,j, the dSBM achieved shorter TTSs than those of a Digital Annealer (DA) (23), which is based on an FPGA implementation of “Digital Annealer’s algorithm” developed from SA and outperformed CPU implementations of SA (25) and parallel tempering (50). (This DA is not the latest version.) Last, for 60-spin and 100-spin problems with all-to-all connectivity, the TTSs of the dSBM are comparable to those of another state-of-the-art machine based on an FPGA implementation of the above-mentioned RBM’s stochastic sampling (28), the size of which, however, is limited to 200 spins. Overall, we conclude that our bSBM and dSBM have achieved remarkably high performance for the present benchmark instances compared to the recently developed machines.

Performance for ultralarge-scale Ising problems

Last, we present the results for two ultralarge-scale Ising problems: a 100,000-spin problem with all-to-all connectivity and a 1,000,000-spin problem with a sparse connectivity of 1% (see section S7 for their detailed definitions for reproduction). Using a GPU cluster with 16 GPUs, we solved these by aSB, bSB, dSB, and our best implementation of SA (see section S10 for details). For comparison, we also solved them by aSB and SA (25) running on a CPU core. Figure 4 shows the results, where the result obtained by the four-GPU implementation of the above-mentioned MA (26) is also shown. The horizontal lines show optimal (dashed) and target (dotted) values estimated using a formula based on statistical mechanics (see section S7 for details) (51, 52).

Fig. 4 Results for ultralarge-scale Ising problems.

(A) Computation times of aSB, bSB, dSB, and our best SA implemented on a 16-GPU cluster, aSB and SA (25) running on a CPU core, and MA implemented on a four-GPU cluster (26) for a 100,000-spin Ising problem with all-to-all connectivity. Each coupling coefficient Ji,j was a randomly chosen 10-bit integer in {−29 + 1,…,29 – 1} (26). Lines and symbols represent average and best values, respectively, for 100 trials (one trial for aSB and SA running on a CPU). (B) Magnification of (A) around the estimated optimal value. (C and D) Corresponding results for a 1,000,000-spin Ising problem with a sparse connectivity of 1%. Each 1% nonzero Ji,j was a randomly chosen 10-bit integer in {−29 + 1,…,29 – 1}.

Figure 4A shows that all three GPU-cluster SBMs outperformed the MA machine (26) in terms of both speed and accuracy. Figure 4 (A and B) also show that the GPU-cluster implementation achieved about 1000 times speedup over a CPU core for aSB but only about 100 times for SA. This difference comes from the higher parallelizability of aSB than that of SA. The GPU-cluster bSBM and dSBM are faster than the GPU-cluster aSBM, because of their algorithmic advantage. Figure 4B shows that the dSBM achieved the closest value to the estimated optimal value, giving it the highest accuracy. As shown in Fig. 4 (C and D), similar results also hold for the 1,000,000-spin Ising problem with sparse connectivity.


In this work, we have proposed two new variants of the SB algorithm, named bSB and dSB, both of which outperform the original aSB in terms of both speed and solution accuracy. dSB allows us to find optimal solutions of large-scale problems, which aSB and bSB cannot achieve. We have implemented 2048-spin-size bSBM and dSBM using single FPGAs. Our benchmarking with TTS and TTT has shown that the bSBM and dSBM can achieve remarkably high performance compared to other recently developed machines. GPU-cluster implementations of bSB and dSB also allow us to find nearly optimal solutions of ultralarge-scale problems with up to one million spins.

Last, we discuss possible future works on SB. First, it is important to check the performance of SB for a broader range of instances than those evaluated in this work, which were chosen to compare our machines with previously developed machines reported in the literature. It is known that a single solver cannot achieve the highest performance for all kinds of instances (53). Thus, we should examine what kinds of instances can be solved well by SB. Second, it is desirable to develop a technique for auto-tuning of parameters in SB, such as a constant c0 and a time step Δt. In this work, we used the definition of c0 based on random matrix theory (30) and also chose the best value of Δt among five values (see Methods for details). Such preliminary search for the best parameter values is often used for benchmarking to evaluate potential performance of solvers. In practical applications, however, it is desirable to determine appropriate parameter values automatically without preliminary search. Last, development of other variants of SB is an interesting possibility. In this work, we have focused on the Ising problem (MAX-CUT). The generalization of SB to broader classes of problem, e.g., combinatorial optimization with higher-order polynomial cost functions, is an interesting next target.


Ballistic simulated bifurcation

In bSB, we numerically solve the Hamiltonian equations of motion given by Eqs. 6 and 7 using the symplectic Euler method (49), as in aSB (30). The updating rule for bSB is as followsyi(tk+1)=yi(tk)+{[a0a(tk)]xi(tk)+c0j=1NJi,jxj(tk)}Δt(15)xi(tk+1)=xi(tk)+a0yi(tk+1)Δt(16)Here, Δt is the time step and tk is the discretized time satisfying t0 = 0 and tk + 1 = tk + Δt. After updating xi, we check whether ∣xi∣ > 1. If so, we replace xi with sgn(xi) and set yi = 0, which implements perfectly inelastic walls at xi = ± 1.

Discrete simulated bifurcation

In dSB, we numerically solve Eqs. 11 and 12. The updating rule for dSB is as followsyi(tk+1)=yi(tk)+{[a0a(tk)]xi(tk)+c0j=1NJi,jsgn[xj(tk)]}Δt(17)xi(tk+1)=xi(tk)+a0yi(tk+1)Δt(18)As in bSB, we replace xi with sgn(xi) and set yi = 0 if ∣xi∣ > 1.

Efficient implementations of bSB and dSB

The bSB and dSB can be implemented more efficiently as follows. [A similar speedup technique has been used for SA (25).]

Storing zi(k)=j=1NJi,jxj(tk), we can rewrite Eq. 15 asyi(tk+1)=yi(tk)+{[a0a(tk)]xi(tk)+c0zi(k)}Δt(19)zi(k)=zi(k1)+j=1NJi,jΔxj(k)(20)where Δxj(k) = xj(tk) − xj(tk − 1). Note that Δxj(k) often becomes zero around the final time. (This is not the case for aSB, making this implementation ineffective.) Hence, the product-sum operation in Eq. 20 can be accelerated by neglecting the terms corresponding to Δxj(k) = 0.

In the case of dSB, instead of Δxj(k), we use Δsj(k) = sgn [xj(tk)] − sgn [xj(tk − 1)] aszi(k)=zi(k1)+j=1NJi,jΔsj(k)(21)Since Δsj(k) becomes zero more often than Δxj(k), this implementation is more effective for speedup in dSB than in bSB.

Parameter setting

In the SB algorithms, the time step Δt and the constants a0 and c0 are appropriately set in advance. In this work, we set the constants as a0 = 1 and c0=0.5JN, where 〈J〉 is defined as J=i,jJi,j2N(N1). This definition of c0 is based on random matrix theory (30). Although this setting of c0 may not be optimal for some instances, this is enough to achieve high performance presented in this work. On the other hand, the setting of Δt is more sensitive to performance. We therefore selected the best value among five values: 0.25, 0.5, 0.75, 1, and 1.25. In Figs. 2A and 4, Δt is set to 0.5 (aSB) or 1 (bSB and dSB). In Fig. 2B, Δt is set to 1.25 for the bSBM and dSBM. The Δt values used in Fig. 3 are provided in table S1. As mentioned in Discussion, automatic setting of c0 and Δt is an important issue but beyond the scope of the present work and so left for future work.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We thank Y. Kaneko, O. Hori, M. Watabiki, M. Iwasaki, M. Ootomo, M. Tomoda, and Y. Izumi for comments and support. Funding: The authors acknowledge no funding in support of this research. Author contributions: H.G. described the dynamical properties of three SBs by producing Fig. 1; estimated optimal solutions for ultralarge instances with T.K.; wrote the manuscript and the Supplementary Materials with data from K.E., Y.H., and R.H.; and supervised the project. K.E. conceived bSB and dSB through discussion with M.S., implemented SBs and SA on a GPU and a 16-GPU cluster, and collected performance data; Y.S. pointed out the mapping of the Ising problem with local fields to one without local fields and convergence to a local minimum of the Ising problem in bSB and dSB; Y.H. and R.H. implemented dSB and bSB, respectively, on an FPGA and collected performance data; M.Y. optimized the FPGA implementations; and K.T. developed the basis of the FPGA implementations, wrote the part of FPGA implementations in the Supplementary Materials, and supervised the FPGA team. Competing interests: H.G., K.E., M.S., Y.S., and K.T. are inventors on a Japanese patent application related to this work filed by the Toshiba Corporation (no. P2019-164742, filed 10 September 2019). The authors declare that they have no other competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested to the authors.

Stay Connected to Science Advances

Navigate This Article