Edge-of-chaos enhanced quantum-inspired algorithm for combinatorial optimization

TL;DR

This work addresses NP-hard Ising-type optimization by extending simulated bifurcation with per-variable nonlinear bifurcation control (GbSB), preserving parallelizable dynamics. GbSB delivers markedly higher solution accuracy, achieving near-100% success on large instances and substantial speedups on FPGA hardware (e.g., 9.6 ms TTS for K2000). The authors uncover an edge-of-chaos mechanism: performance peaks near the boundary between regular and chaotic dynamics, enabling chaos-assisted search. These findings, along with a highly parallel GbSB FPGA machine, suggest physics-inspired dynamical systems can tackle intractable combinatorial problems with unprecedented speed, albeit with caveats and potential need for hybrid approaches for certain instances.

Abstract

Nonlinear dynamical systems with continuous variables can be used for solving combinatorial optimization problems with discrete variables. Numerical simulations of them are also useful as heuristic algorithms with a desirable property, namely, parallelizability, which allows us to execute them in a massively parallel manner, leading to ultrafast performance. However, the dynamical-system approaches with continuous variables are usually less accurate than conventional approaches with discrete variables such as simulated annealing. To improve the solution accuracy of a quantum-inspired algorithm called simulated bifurcation (SB), which was found from classical simulation of a quantum nonlinear oscillator network exhibiting quantum bifurcation, here we generalize it by introducing nonlinear control of individual bifurcation parameters and show that the generalized SB (GSB) can achieve surprisingly high performance, namely, almost 100% success probabilities for some large-scale problems. As a result, the time to solution for a 2,000-variable problem is shortened to 10 ms by a GSB-based machine, which is two orders of magnitude shorter than the best known value, 1.3 s, previously obtained by an SB-based machine. To examine the reason for the ultrahigh performance, we investigated chaos in the GSB changing the nonlinear-control strength and found that the dramatic increase of success probabilities happens near the edge of chaos. That is, the GSB can find a solution with high probability by harnessing the edge of chaos. This finding suggests that dynamical-system approaches to combinatorial optimization will be enhanced by harnessing the edge of chaos, opening a broad possibility to tackle intractable combinatorial optimization problems by physics-inspired approaches.

Figures (5)

• Figure 1: Numerical results of the GbSB. (a) Success probability, $P_\mathrm{S}$, of finding the best known cut value of K$_{2000}$ by the GbSB with ${\Delta_t = 1.25}$ and various values of $M$ and $A$. See Appendix \ref{['Appendix1']} for the setting of $c$. The lower figure is the enlargement of the part enclosed by the dashed rectangle in the upper figure. $P_\mathrm{S}$ was evaluated by 1,000 times repetition with different initial conditions. (b) and (c) Results corresponding to (a) for an instance of the 700-spin Ising problem and G6, respectively. The time step is set as ${\Delta_t = 1.25}$ and ${D_t = 1.25}$ [see Eq. (\ref{['eq13']}) in Appendix \ref{['Appendix1']}], respectively. See Appendix \ref{['Appendix1']} for the setting of $c$. (d) $P_\mathrm{S}$ in a (upper figure) and normalized distance, $\delta (t_M)$, at the final time (lower figure) for K$_{2000}$ when ${M=21,500}$ in (a). $\delta (t_M)$ was evaluated by averaging 100 results with different initial conditions. The regions highlighted in blue show the condition that $P_\mathrm{S}$ becomes particularly high. The dotted line in the lower figure shows $\delta (t_M )=1/\sqrt{2}$, which indicates chaos. (e) and (f) Results corresponding to (d) for the 700-spin instance when ${M=3,000}$ in (b) and G6 when ${M=1,000}$ in (c), respectively. (g--i) $P_\mathrm{S}$ for K$_{2000}$, the 700-spin instance, and G6, respectively, for different values of the time step $\Delta_\mathrm{t}$ and different final times $t_M=\Delta_\mathrm{t} M$. $A$ was set to 0.2, 0.25, and 0.2, respectively. (j--l) $\delta (t_M)$ corresponding to (g--i), respectively.
• Figure 2: GbSB-based machine with an FPGA. (a) Block diagram of the GbSB-based machine, which has a circulative structure corresponding to iteration of updating position and momentum: $y_i (t_{m+1} )=y_i (t_m )+\delta y_i (t_m ) + \Delta y_i (t_m )$ and $x_i (t_{m+1} )=x_i (t_m )+\delta x_i (t_m )$, where $\Delta y_i (t_m ) = \mathrm{JX}[\mathbf{J}_i, \mathbf{x}(t_m)] = c \Delta_\mathrm{t} \sum_{j=1}^N J_{i,j} x_j (t_m)$, $p_i (t_{m+1} )=\mathrm{FP}[p_i (t_m ), x_i (t_m )] = p_i (t_m ) - [1-A x_i (t_m)^2 ] p_i (t_m )/(M-m)$, $\delta y_i (t_m ) = \mathrm{FX}[p_i (t_{m+1} ), x_i (t_m )] = -p_i (t_{m+1} ) x_i (t_m ) \Delta_\mathrm{t}$, and $\delta x_i (t_m ) = \mathrm{FY}[y_i (t_{m+1} )] = y_i (t_{m+1} ) \Delta_\mathrm{t}$. The JX, FP, FX, FY modules in the block diagram correspond to the JX, FP, FX, FY functions, and the wall module executes the following conditional operation: $(x_i (t_{m+1} ),y_i (t_{m+1} )) \leftarrow \left\{ (\mathrm{sgn} (x_i (t_{m+1} )), 0),\mathrm{if}~|x_i (t_{m+1} )|>1(x_i (t_{m+1} ), y_i (t_{m+1} )),\mathrm{otherwise} \right.$. The computational precision of the JX module is 16-bit fixed point, while that for the remaining modules is 32-bit floating point. (b) Block diagram of a many-body interaction (JX) module. The parallelization parameters ($P_\mathrm{r}$, $P_\mathrm{c}$, and $P_\mathrm{b}$) are illustrated in a and b. (c) Layout of the circuit modules in the FPGA, where the routing congestion is shown as a heatmap with the regions for the JX, PXY, and BSP (board support package) modules being indicated by dashed lines.
• Figure 3: Times to solution (TTSs) with a GbSB-based FPGA machine. The second row (${N=700}$) is the median of the TTSs for 100 random instances of the 700-spin Ising problem. See Appendix \ref{['Appendix1']} for the detailed parameter settings and the GbSB-based machine (GbSBM). The results with the dSB-based machine (dSBM) are cited from ref. 21.
• Figure 4: Average cut value $C_{\mathrm{ave}}$ normalized by the best know value $C_{\mathrm{best}}$ corresponding to Figs. \ref{['fig1']}(a--c), respectively.
• Figure 5: Parallel processing of the computation in a GbSB time-evolution step. (a) Timing chart of the operation of JX and PXY modules. (b) A schematic showing how the many-body interaction computations, $\sum_{j=1}^N J_{i,j} x_j (t_m)$, are processed in parallel.