Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems

SB algorithm

The N-spin Ising problem without external magnetic fields is formulated as follows. A dimensionless Ising energy is given byEIsing(s)=−12∑i=1N∑j=1NJi,jsisj(1)where s_i denotes the ith Ising spin, which takes a value of 1 or −1, s = (s₁ s₂ … s_N) is the vector representation of a spin configuration, and J_i,j(= J_j,i) is the coupling coefficient between the ith and jth spins (J_i,i = 0). The problem is to find a spin configuration minimizing the Ising energy. This problem is mathematically equivalent to a famous combinatorial optimization problem named MAX-CUT (5, 17, 18, 39): divide the nodes of a weighted graph into two groups maximizing the total weight (called “cut value”) of the edges cut by the division. By setting the coupling coefficients as J_i,j = −w_i,j (w_i,j = w_j,i is the weight between the ith and jth nodes), we can solve MAX-CUT using Ising machines (see also the Supplementary Materials) (5, 17, 18). Since the total number of spin configurations is 2^N, it is extremely difficult to find exact solutions for large N. Unless the topology of the graph has some special structure, such as two-dimensional lattices, this problem is known to be nondeterministic polynomial time (NP)–hard (2, 3). Hence, it is strongly believed that there exists no efficient exact algorithm for this problem. In practice, one often uses heuristic algorithms, such as SA, to find approximate but practically useful solutions as fast as possible. Here, we propose SB as a new option for these heuristic algorithms. [For obtaining exact solutions of small-size problems, the SA machine called “Digital Annealer” in (29) may be the fastest so far.]

For the Ising problem, a new kind of quantum adiabatic optimization using a network of Kerr-nonlinear parametric oscillators (KPOs) has recently been proposed (30–34). The quantum mechanical Hamiltonian in this approach is given by (30, 34)Hq(t)=ℏ∑i=1N[K2ai†2ai2−p(t)2(ai†2+ai2)+Δiai†ai]−ℏξ0∑i=1N∑j=1NJi,jai†aj(2)where ℏ is the reduced Planck constant, ai† and a_i are the creation and annihilation operators, respectively, for the ith oscillator, K is the positive Kerr coefficient, p(t) is the time-dependent parametric two-photon pumping amplitude, Δ_i is the positive detuning frequency between the resonance frequency of the ith oscillator and half the pumping frequency, and ξ₀ is a positive constant with the dimension of frequency. The initial state of each KPO is set to the vacuum state, and the pumping amplitude p(t) is gradually increased from zero to a sufficiently large value. The constant ξ₀ is set to a sufficiently small value such that the vacuum state is the ground state of the initial Hamiltonian (30, 34). Then, each KPO finally becomes a coherent state with a positive or negative amplitude via a quantum adiabatic bifurcation (30, 34), and the sign of the final amplitude for the ith KPO provides the ith Ising spin of the ground state of the Ising model, which is guaranteed by the quantum adiabatic theorem (30, 34).

The corresponding classical Hamiltonian system is derived from a classical approximation where the expectation value of a_i is approximated as a complex amplitude x_i + iy_i (where i denotes the imaginary unit) (30, 34). Here, the real and imaginary parts, x_i and y_i, are a pair of canonical conjugate variables, such as a position and a momentum, for the ith KPO. The classical mechanical Hamiltonian and the equations of motion for this classical system are given by (30, 34)Hc(x,y,t)=∑i=1N[K4(xi2+yi2)2−p(t)2(xi2−yi2)+Δi2(xi2+yi2)]−ξ02∑i=1N∑j=1NJi,j(xixj+yiyj)(3)x.i=∂Hc∂yi=[K(xi2+yi2)+p(t)+Δi]yi−ξ0∑j=1NJi,jyj(4)y.i=−∂Hc∂xi=−[K(xi2+yi2)−p(t)+Δi]xi+ξ0∑j=1NJi,jxj(5)where the dots denote differentiation with respect to time t. It was empirically found that this classical system can also find good approximate solutions with considerably high probability, where the sign of the final x_i provides the ith spin of the solutions (30, 34). This result suggests a new approach to the Ising problem. Unlike quantum computers, we do not need to build a real machine described by the above Hamiltonian, because, instead, we can simulate such a machine efficiently using classical computers.

However, Eqs. 4 and 5 are not suitable for fast numerical simulation. By disregarding some terms proportional to the momenta y, which varies around zero (30, 34), we simplify Eqs. 4 and 5 as followsHSB(x,y,t)=∑i=1NΔ2yi2+V(x,t)=∑i=1NΔ2yi2+∑i=1N[K4xi4+Δ−p(t)2xi2]−ξ02∑i=1N∑j=1NJi,jxixj(6)x.i=∂HSB∂yi=Δyi(7)y.i=−∂HSB∂xi=−[Kxi2−p(t)+Δ]xi+ξ0∑j=1NJi,jxj(8)where V(x, t) denotes the potential energy, and all the detunings have been assumed to be the same value Δ. The resultant classical system is a network of Duffing oscillators (35, 40) with mass Δ⁻¹ and coupling coefficients {ξ₀J_i,j}. [Similar simplification for the CIM simulation has recently been reported in (41).]

Hereafter, we treat all the parameters with the dimension of frequency (correspondingly, also time) as dimensionless quantities because we use Eqs. 7 and 8 as an algorithm for solving mathematical problems, and therefore, the physical meanings of the parameters are unimportant. x and y are also dimensionless by definition.

The SB algorithm is based on Eqs. 7 and 8. Note that the new Hamiltonian in Eq. 6 is separable with respect to positions and momenta, unlike the previous one in Eq. 3. Therefore, Eqs. 7 and 8 can be solved by the explicit symplectic Euler method (see also Methods) (36) instead of implicit methods or, e.g., the explicit fourth-order Runge-Kutta method. The simplicity of the explicit symplectic Euler method is crucial for hardwiring of the SB algorithm with custom circuits, such as FPGAs. Its stability is also important, because this allows one to set the time step to a large value, which results in high-speed simulation. [The symplectic Euler method is first-order accurate. Higher-order explicit symplectic methods, such as the Störmer-Verlet method (36), can also be used, but this is not effective not only for the speed and stability of the present simulation but also for the solution accuracy for the Ising problem. Hence, we adopt the simplest first-order method. See Methods for details.]

The SB algorithm is as follows. All the variables, x and y, are initially set around zero. Gradually increasing p(t) from zero, we solve Eqs. 7 and 8 numerically by the explicit symplectic Euler method. (For further speedup, we modify the method from its standard form. See Methods for details.) The sign of the final x_i provides the ith spin of an approximate solution of the Ising problem.

SB has further algorithmic advantages: Eqs. 7 and 8 have only one product-sum operation with respect to the coupling matrix J, which is the most computation-intensive part in this algorithm, unlike Eqs. 4 and 5 including two such product-sum operations; we can update the variables simultaneously, and therefore, we can fully exploit massively parallel processing for accelerating SB; and SB includes only addition and multiplication, and no probabilistic processes, and therefore can be easily hardwired, e.g., with FPGAs.

Mechanism of SB

The mechanism of SB is qualitatively explained as follows. When p(t) is gradually increased from zero to a sufficiently large value p_f, each oscillator exhibits a bifurcation with two stable branches whose positions are approximately given by si(pf−Δ)/K (s_i = ±1) (34). Correspondingly, the potential energy finally has 2^N minima corresponding to 2^N spin configurations, which are approximately given byV(s(pf−Δ)/K,t)=−N(pf−Δ)24K−ξ02pf−ΔK∑i=1N∑j=1NJi,jsisj(9)

Note that the first term is independent of s, and the second term is proportional to the Ising energy regarding s as a spin configuration. In the adiabatic evolution with slowly varying p(t), the state will follow one of the potential energy minima near to the present state. Assuming that minima with relatively low final energies appear earlier and are lower during the adiabatic evolution than those with relatively high final energies, we expect that the state will arrive at a minimum with a relatively low final energy. Thus, we will find a spin configuration with a low Ising energy as an approximate solution.

Instead of showing the general validity of the above qualitative description, here, we show the situation of a simple example: the two-spin Ising model with ferromagnetic coupling (J_1,2 = J_2,1 = 1), the ground states of which are (up, up) and (down, down). This will give us an intuitive understanding of SB, which is enough to use SB as a heuristic algorithm.

Figure 1 (A and B) shows the simulation results for this problem with a linearly increased pumping amplitude p(t) = ε_pt. As shown in this figure, after two bifurcations, the potential energy finally has four minima corresponding to the four configurations of the two-spin Ising model.

• Download high-res image
• Open in new tab
• Download Powerpoint

Before the first bifurcation, the potential energy has a single minimum at the origin, around which the trajectory circulates, as shown in the top panel of Fig. 1B. After the first bifurcation, two stable fixed points (potential energy minima) corresponding to the ground states, (up, up) and (down, down), appear earlier than the other two minima corresponding to (up, down) and (down, up). The state adiabatically follows one of the two minima that appear earlier. Thus, SB successfully finds a ground state (down, down).

This adiabatic evolution is explained as follows. First, the total energy E(t) and the potential energies at the minima, V_±,±(t) (+ and − correspond to up and down, respectively), are formulated from their time derivatives asE(t)=E0−εp2∫0t|x(t′)|2dt′(10)V±,±(t)=−εp2∫0t|X±,±(t′)|2dt′(11)where E₀ is the initial total energy, which is assumed to be small, and X_±,±(t) denote the positions of the minima. Before the corresponding bifurcation, X_±,±(t) are defined as zero. Note that the total energy monotonically decreases because of the second term in Eq. 10 describing the negative work done by the time-dependent pumping term. After the first bifurcation, X_+,+(t) and X_−,−(t) leave the origin, and x(t) starts to circulate around them, instead of the origin. This leads to the increase of |x(t)| and therefore the further decrease of the total energy. As a result, the energetically allowable region defined by V(x, t) − E(t) ≤ 0 (inside the loops in Fig. 1B) separates into two parts, as shown in the middle panel of Fig. 1B. In the present simulation, x(t) is confined in the allowable region around X_−,−(t). After that, x(t) circulates around X_−,−(t), which leads to ∫|x(t)|²dt ≈ ∫|X_−,−(t)|² dt. From this and Eqs. 10 and 11, the difference between E(t) and V_−,−(t) is kept constant at a small value. Thus, x(t) is kept in the allowable region around X_−,−(t), as shown in the bottom panel of Fig. 1B. In summary, SB finds a ground state by adiabatically following one of the stable fixed points (potential energy minima) corresponding to the ground states that appear at the first bifurcation point.

The above mechanism of SB suggests that solutions obtained by SB are determined at the first bifurcation point. Since x(t) at the first bifurcation point is given by the maximum-eigenvalue eigenvector of the coupling matrix J (see the Supplementary Materials), SB may provide the signs of the elements of the eigenvector as an approximate solution, which is actually an approximate solution obtained by a continuous reduction method (30, 34). As shown in the next subsection, however, SB can provide much better solutions than this. This fact implies that there may exist another mechanism of SB; that is, SB may continue to search better solutions after the first bifurcation.

To investigate this new mechanism of SB, we also did another simulation where p(t) is given in the top panel of Fig. 1C. The results are shown in Fig. 1 (C and D). In this case, the initial value of p(t) is higher than the two bifurcation points, and therefore, there are four potential energy minima even at the initial time. While p(t) is constant, x(t) moves in the energetically allowable region including the four minima, as shown in the top panel of Fig. 1D. After that, the allowable region separates into four parts, as shown in the middle panel of Fig. 1D. Then, x(t) is confined in one of the four regions. Eventually, SB find a ground state (up, up) successfully.

Since x(t) also moves around the other potential energy minima, as shown in the top and middle panels of Fig. 1D, the failure probability may be finite. To evaluate the success probability, we repeat this simulation 10⁴ times with different initial conditions, where x₁(0) = x₂(0) = 0 and y₁(0) and y₂(0) are set randomly from the interval (−0.1, 0.1). As a result, we obtain ground states with a probability of about 0.83. This result indicates that lower energy states (better solutions) are obtained with higher probabilities even in the initial condition where the energetically allowable region includes multiple potential energy minima.

This interesting result can be explained by assuming the ergodicity of this nonlinear Hamiltonian system, implying that a trajectory will eventually visit all states in the phase-space allowable region (36). This is suggested by the top panel of Fig. 1D, where x(t) moves around two low minima more frequently than two high minima. From the ergodicity, the success probability may be proportional to the phase-space volume of the allowable region satisfying x₁x₂ > 0. The success probability numerically estimated by the phase-space volume for the middle panel of Fig. 1D is about 0.841, which is remarkably close to the actual success probability of 0.83. This numerical fact supports the mechanism based on the ergodicity. In general, the ergodicity suggests that x(t) moves around low potential energy minima in longer time than high minima, which results in lower Ising energies. In the next subsection, we provide additional numerical evidence for this mechanism; that is, SB can find much better solutions for a 2000-spin Ising problem than that obtained at the first bifurcation. Thus, we conclude that SB finds lower minima among many minima harnessing the ergodicity. [Very recently, a new approach to the Ising problem using a classical Hamiltonian system with “classical spins” mimicking quantum annealing has been proposed (42), where a fixed point is tracked by the technique called shortcut to adiabaticity. Although this approach uses classical Hamiltonian dynamics as well, it is quite different from SB, because it will not exploit the ergodic search but just track a fixed point. In addition, it has been unclear so far whether it can achieve a high-speed computation with its complicated equations.]

Solving an all-to-all connected 2000-node MAX-CUT problem by a single-FPGA SB machine

From the advantages and properties of SB, we expect that SB will be useful for approximately solving large-scale Ising or MAX-CUT problems with dense connectivity as fast as possible. To check this expectation, we compare our SB machine with a state-of-the-art CIM (17, 19), because this machine has achieved outstanding performance for these problems. For this comparison, we solved an all-to-all connected 2000-node MAX-CUT problem named K₂₀₀₀, which was solved by the CIM (17, 19).

We implemented an SB-based 2000-spin Ising machine with a single FPGA. The product-sum operation, ∑j=1NJi,jxj, in Eq. 8, which is the most computation-intensive part in SB, is performed efficiently by massively parallel processing. J_i,jx_j terms of up to 8192 are processed at a single clock, which is about four times more than the total spins. See the Supplementary Materials for details.

One of the main results in (17) is that the CIM reached a target energy given by the Goemans-Williamson semidefinite programming (GW-SDP) algorithm (17, 38) in about 70 μs in the best case among 100 trials, which is about 50 times shorter than that by the SA highly tuned in (17). [The SA in (17, 19) is very fast by putting all the elements of the coupling matrix J on a cache of a single core. Thus, parallel computing with multiple cores cannot accelerate the SA because of communication overheads.] In comparison with this, our single-FPGA SB machine reaches the GW-SDP value in only 58 μs even in the worst case among 100 trials, as shown in Fig. 2A.

• Download high-res image
• Open in new tab
• Download Powerpoint

Another main result in (17) is that the CIM took only 5 ms to obtain as good solutions (large cut values) as the highly tuned SA obtained in 50 ms, where the cut value of the MAX-CUT problem is given by −EIsing2−14∑i=1N∑j=1NJi,j with Ising parameters (see the Supplementary Materials). Figure 2B shows that our SB machine takes only 0.5 ms to obtain better solutions on average than both the CIM in 5 ms and the SA in 50 ms.

In Fig. 2A, we also show the energy obtained by the Hopfield neural network (HNN) of two-state neurons (19, 43). The HNN is the most naïve local search algorithm, and therefore, more sophisticated algorithms should outperform the HNN. (The HNN is equivalent to SA at zero temperature; see the Supplementary Materials.) Note that the HNN energy is lower than the GW-SDP energy. Our SB machine reaches the HNN energy about 13 times faster than the CIM on average. From these results, we conclude that our single-FPGA SB machine is about 10 times faster than the CIM developed in (17).

As mentioned in the last subsection, at the first bifurcation point, SB may find an approximate solution determined by the maximum-eigenvalue eigenvector of the coupling matrix J (30, 34). In the case of K₂₀₀₀, this approximate solution is E_Ising = −55724 or the cut value of 27,342. This is much worse than those actually obtained by SB, as shown in Fig. 2. This fact supports the mechanism of SB based on the ergodicity discussed in the last subsection.

Solving an all-to-all connected 100,000-node MAX-CUT problem by a GPU-cluster SB machine

We have further solved an all-to-all connected 100,000-node MAX-CUT problem with continuous weights set randomly from the interval [−1, 1], which is named M₁₀₀₀₀₀. [See the Supplementary Materials for its definition with a pseudo-random number generator in (44).] Note that in M₁₀₀₀₀₀, there are about 5 billion edges, the data size of which is too large for FPGAs. We instead use two large computational systems and compare the performance from each. System 1 is a GPU cluster using eight GPUs communicating with one another via a fast link. System 2 is a PC cluster composed of multiple nodes with two 14-core processors each communicating with one another via a fast link. We also solved M₁₀₀₀₀₀ by the HNN and the SA developed by us based on (22). Our SA is remarkably fast, e.g., about 100 times faster than the software provided by (22). See the Supplementary Materials for details.

The comparison of computation times for SB and SA is shown in Fig. 3A. The fastest implementation is the GPU cluster for SB and the 50-core PC cluster for SA. This difference comes from the fact that the massively parallel processing of the GPU cluster is effective for SB, but not for SA, due to the lack of parallelizability of the spin update and resultant large communication overheads. Figure 3B shows the computation times of SB and SA in their fastest implementations. The dotted lines roughly give the envelopes of the time-evolution curves. This means that the dotted lines indicate the best cases with respect to N_step and N_sweep. From this comparison at the same Ising energies below the HNN value, it turns out that the GPU-cluster SB machine is about 10 times faster than our fastest SA or 1000 times faster than the software provided by (22). This demonstrates that SB allows one to accelerate large-scale combinatorial optimization without specially designed hardware devices or custom circuits.

• Download high-res image
• Open in new tab
• Download Powerpoint