Bifurcation behaviors shape how continuous physical dynamics solves discrete Ising optimization

Simulating physical dynamics to solve hard combinatorial optimization has proven effective for medium- to large-scale problems. The dynamics of such systems is continuous, with no guarantee of finding optimal solutions of the original discrete problem. We investigate the open question of when simulated physical solvers solve discrete optimizations correctly, with a focus on coherent Ising machines (CIMs). Having established the existence of an exact mapping between CIM dynamics and discrete Ising optimization, we report two fundamentally distinct bifurcation behaviors of the Ising dynamics at the first bifurcation point: either all nodal states simultaneously deviate from zero (synchronized bifurcation) or undergo a cascade of such deviations (retarded bifurcation). For synchronized bifurcation, we prove that when the nodal states are uniformly bounded away from the origin, they contain sufficient information for exactly solving the Ising problem. When the exact mapping conditions are violated, subsequent bifurcations become necessary and often cause slow convergence. Inspired by those findings, we devise a trapping-and-correction (TAC) technique to accelerate dynamics-based Ising solvers, including CIMs and simulated bifurcation. TAC takes advantage of early bifurcated “trapped nodes” which maintain their sign throughout the Ising dynamics to reduce computation time effectively. Using problem instances from open benchmark and random Ising models, we validate the superior convergence and accuracy of TAC.

Proof. Consider a set of compact subsets of Euclidean space R 2n : R = {d(s 1 , s 2 , . . . , s n ) = I s1 × I s2 × · · · × I sn × I n a | s i ∈ {−1, 0, +1}} , Daniel Ebler: ebler.daniel@huawei.com Jie Sun: j.sun@huawei.com where I a is defined as Then for d(s 1 , s 2 , . . . , s n ) ∈ R, consider a mapping F : d(s 1 , s 2 , . . . , s n ) → R 2n such that F (x, y) = (F 1 (x, y), F 2 (x, y), . . . , F 2n (x, y)) where x ∈ I s1 × I s2 × · · · × I sn and y ∈ I n a , and where g k (x) = x 3 −(−1+p)x and x ∈ I s k . Note that it can be easily verified that ∀s k ∈ {−1, 0, +1}, we have g k is monotonic in I s k and g k : I s k → I a is bijective, so its inverse g −1 k : I a → I s k exists and is continuous. Then we can conclude that F is a continuous mapping. Next we will show that if p > 1 + 1 For (1), it is valid by definition, as I si ∩ I s i = ∅, for s i = s i . For (2), consider (s 1 , s 2 , . . . , s n ) ∈ {−1, 0, +1} n , and let x ∈ I s1 × I s2 × · · · × I sn and y ∈ I n a . For k = 1, 2, . . . , n, we have For k = n + 1, n + 2, . . . , 2n, if p > 1 + 1 and hence F k (x, y) ∈ I a , for k = n + 1, n + 2, . . . , 2n. Thus we have shown that if p > 1 +  (1) and (2) and Brouwer's Fixed-Point Theorem [1], for any d(s 1 , s 2 , . . . , s n ) ∈ R, there exists a fixed point (x , y ) ∈ d(s 1 , s 2 , . . . , s n ) such that F (x , y ) = (x , y ), or equivalently, . . . , n. As |R| = 3 n by definition, we conclude that there are at least 3 n solutions of CIM dynamics Supplementary Equation (2). By the Finiteness Theorem (or specifically, by checking the condition on Gröner basis of the system of fixed points equation of Supplementary Equation (2)), this solution set of the system is zero-dimensional (as for each i, the Gröner basis has an element whose leading monomial is a power of x i ), hence by Bézout's Theorem [2], there are at most 3 n solutions. Therefore, we can further conclude that for any (s 1 , s 2 , . . . , s n ) ∈ {−1, 0, +1} n , there exists exactly one x ∈ I s1 × I s2 × · · · × I sn such that where D = max l j |ξG lj |, then there exists a stable fixed point of CIM dynamics (2) whose binarization corresponds to the ground states of Ising problem Supplementary Equation (1).
Next, we show that x is stable. The Jacobian matrix of the dynamics Supplementary Equation (2) at x is given as By Courant-Fischer min-max Theorem and Lemma 1, its maximum eigenvalue λ max (J| x ) satisfies where λ max (·) is the maximum eigenvalue of the given matrix. Notice that and therefore we have showing that x is a stable fixed point.

Supplementary Note 2 -Emergence of non-zero stable state at the first bifurcation
In this section, we derive the state of Coherent Ising Machine (CIM) after the first bifurcation point p = p 0 explicitly and prove its stability. The dynamical equations of the n-dimensional system x = (x 1 , . . . , x n ) T are given as and the corresponding potential function is For control parameter p < p 0 the variables assume the trivial state x 0 = 0. When this state gets unstable, the system bifurcates into two separate branches, corresponding to the newly arising fixed points in Supplementary Equation (20). Depending on the problem, one or both of the branches can be stable, and generally (if no particular symmetries or structures are present) lead to different energy levels.
The initial state x 0 gets unstable when the Jacobian J| x = −∇ 2 U (x) has positive maximum eigenvalue λ max (J| x ). Explicitly, the Jacobian reads and yields for the trivial state Analyzing the spectrum of J| x0 , we find the parameter p 0 as the minimum p satisfying Consequently, the trivial solution becomes unstable at p 0 = 1−ξλ max (G) and the state x 0 bifurcates into x 1 ∝ v max (G) (see main text). Next, we show that x 1 is a stable state. Stability at a point x is equivalent to the Jacobian J| x = −∇ 2 U (x) having positive maximum eigenvalue λ max (J| x ). To perturb J| vmax around the bifurcation point p 0 , we define J| x1 = J| x0 + δJ, with δJ = −3diag(δ 2 1 , ..., δ 2 n ) (see Supplementary Equation (22)). Then, we find [3] following from λ max (J| x0 ) = 0 and negative definiteness of δJ. Hence, x 1 corresponds to a stable local minimum. In the following, we further derive the evolution of amplitudes of fixed points when p is close to p 0 . Given that the bifurcation direction of the CIM is proportional to the maximum eigenvector v max of the interaction matrix G, the stable CIM state at pump rate close to bifurcation point p = p 0 can be approximated by x = av max . Then the parameter dependence of the amplitude of the fixed point ||x|| 2 can be obtained as where we denote K = I − 1 λmax(G) G and D = diag(v max,1 , . . . , v max,n ). Note that K is a singular matrix with kernel spanned by v max . Here without loss of generality, we assume the dimension of the kernel is one. As K is a symmetric matrix, its singular value decomposition can be written as where e 1 = (1, 0, 0, . . . , 0) T . Notice that the (i, j)-th element of an inverse matrix M ij is given as of M by removing j-th row and i-th column. As the pump rate is close to bifurcation threshold, we have a → 0, and cof(S + 3a 2 U T D 2 V, i, j) = cof(S, i, j), and therefore, Notice that so Supplementary Equation (29) can be simplified as By plugging Supplementary Equation (31) into Supplementary Equation (27), we have

Supplementary Note 3 -Bifurcation analysis of the embedded potential
In this section we analyze the dynamics of the local minima points of the potential U (x) during the bifurcation cascade of the system. This complements the discussion on the bifurcation behavior of the fixed points in the main text. Supplementary Figure 1 shows the dynamics of a local minimum under increasing pump rate p for an outlier with n = 5 variables. Until the first bifurcation threshold p 0 the value of the minima is constant at U (x 0 ). When reaching p 0 (solid grey line), some of the nodes bifurcate -leading to new emerging minima in U -while others stay exactly zero. Hence, a mapping to the landscape of the discrete optimization problem is not possible at this point. Note that this is generally not the only possible scenario: all variables can be -close to zero (swing nodes) after the bifurcation, but not exactly zero. In such a setting, the nodes have bifurcated but their magnitudes remain less than -which enables a mapping onto the QUBO landscape.
At Supplementary Figure 1: Evolution of the potential U of local minima points after the bifurcation in an Ising model with n = 5 variables. In this example, the decision point p * > p0 is larger than the bifurcation threshold, causing the problem to be an outlier.
Supplementary Figure 3 illustrates the same analysis for a larger system (n = 20). For the sake of visibility, only the evolution of the lowest minima of U at each value p is depicted. This example is slightly more involved. Different from the previous case of n = 5, the first bifurcation at p = p 0 causes all nodes to drift away from zero as shown in Supplementary Figure 3(b). In this instance, no node stays exactly zero. Consequently, the mapping of the minima of U (x) between p 0 and p * onto the real problem is possible -in contrast to the case of n = 5 where some variables remained exactly zero. Nevertheless, variables x i with values close to zero indicate that they are swing nodes, since their small values result from cancellation of opposite fields of their neighbors. The binarization of their values leads to random assignments. This gives, ultimately, a mapping of U (x) onto local minima of H Ising -hence, the green solid line in Supplementary Figure 3. For this example, the middle orange line (the most centric line around zero) represents state of a node whose value remains -close to zero -no escape from the -neighborhood around zero after p 0 , but it is not until the pump rate reaches p * (marked by the vertical grey dashed line) that it escapes and crosses the zero line (we call it a crossover) and flip its sign from −1 to +1 leading to the exact

mapping.
Until the point p = p * , several further bifurcations and crossovers trap some of the nodes close to zero. For the potential U (x), this means that new minima (stable and metastable) emerge, while existing ones evolve in the landscape (see Supplementary Figure 4). After the point p = p * is reached, all nodes are trapped or otherwise finished their crossover. One minimum in U (solid blue line) evolves, with the corresponding state x * mapped to the ground state of the the discrete problem. Then based on Theorem 1, we can further show that if all nodes of a system exhibit a synchronized bifurcation -meaning that the magnitudes |x 1,i |, ∀i = 1, 2, ..., n of all nodes are sufficiently large at the first bifurcation point p 0 , then the mapping between the continuous dynamics and the original optimization problem is exact.

Corollary 1. Suppose that the conditions stated in Theorem 1 hold, it follows that if
then σ 1 = sign(x 1 ) ∈ {−1, +1} n is a ground state of the Ising problem.

Supplementary Note 5 -Retarded bifurcation detection
In this section we briefly describe the observations of retarded bifurcations. CIMs are commonly capable to decide a certain number of variables efficiently (at p = p 0 ). For such variables i, the magnitude |x i | is significantly larger than zero. Consequently, further evolution of the system (across additional bifurcation points) are not likely to change the sign of the variables anymore.
Variables j for which |x j | < , with 1 can require longer (up to exponential) run-times to be trapped. These variables can be the cause of higher computational complexity of the respective optimization tasks. Hence, such retarded bifurcations are usually the cause for CIMs to output sub-optimal solutions, as the runs often do not reach p = p * or the CIM fails to stay at the state with the lowest potential during the evolution.
We exhaustively investigate the maximum eigenvectors of all possible instances with n = 16 of Ising problems on cubic graphs -i.e. graphs for which each node shares three links to other nodes (this set of graphs was previously studied, for instance, in Ref. [4]). The connections are G ij = −1 so the Ising problems correspond to the max-cut problems on the cubic graphs. Supplementary Figure 5 shows the normalized counts of the magnitudes |v max,i | for all 4060 possible problems. It can be seen from Supplementary Figure 5(a) that for all problems with p * = p 0 , i.e., σ = sign(v max ) is the ground state, the detected number of nodes close to zero is low. In contrast, problems subject to retarded bifurcations (with p * > p 0 ) exhibit a proportionately dominant number of nodes close to zero -see Supplementary Figure 5 Supplementary Note 6 -Swing nodes and correctness of the solution at first bifurcation p 0 In this section we illustrate the statistics of the swing nodes and the degree of synchronization as a function of the size of graphs n in Supplementary Figure 6. It can be seen that the presence of retarded bifurcation (low degree of synchronization) and swing nodes become more significant when n increases, and they lead to increasing deviations (decreasing correctness) of the binary solution corresponding to the first non-trivial stable state of CIM from the correct solution (see Supplementary Figure 7).
Supplementary Figure 6: Statistics of the swing nodes and the degree of synchronization for graphs with different number of nodes n. For each n, 100 fully connected graphs with random binary weights (±1) are generated to obtain the statistics. (a) The percentage of graph instances that contain at least one swing node at first bifurcation p = p0 (blue bars) and the average percentage of swing nodes (orange curve) as a function of graph size n. Here the criterion for classifying nodes by first computing the state at the first bifurcation point as x1 ∝ vmax, and then defining swing nodes as those j-th element such that |x1,j| < , where the threshold is = 0.5||x1||/ √ n. This figure shows that while the absolute number of swing nodes tend to increase with graph size, their percentage tends to saturate thus providing an opportunity to utilize the remainder (trapping nodes) when searching for an optimal solution. (b) Average degree of synchronization α 2 (x1) as a function of n, where the mean (solid curve) and one standard deviation (shaded area) are both shown. The saturation α 2 (x1) corresponds to the saturation of swing nodes fraction as in (a).
Supplementary Figure 7: The correct fraction r using the binarized state at the first bifurcation x1 ∝ vmax as the function of the graph size n. Let s = sign(x1) = sign(vmax) and denote the ground state configuration as σ * , the correct fraction is the fraction of binarized nodes agreeing with the ground state configuration, taken to be either σ * or −σ * due to the inversion symmetry of the objective function. It is computed as r = max (|{i : si = σ * i }|, |{i : si = −σ * i }|), where | · | represents cardinality of the set. For graphs with size n > 10, we estimate σ * by running SB, bSB, as well as dSB each with 100 random initializations for an extended time and choosing the binary configuration that gave the lowest Ising energy. Average r(n) is shown as the solid curve and shaded by one standard deviation, with a logarithmic fit denoted using the dashed line. In this figure, for each graph size n, a total of 20 fully connected graphs with random binary weights (±1) are generated to obtain the statistics.

Supplementary Note 7 -Trapping-and-correction (TAC) approach
In this section we define the trapping-and-correction (TAC) approach proposed in the main text and elaborate on the single steps in detail.
The approach takes a system state and a cutoff parameter = k||x||/ √ n as inputs. The integer k is fixed by the experimenter. Suppose the system dynamics stops at an instantp and given the state xp of the CIM or other dynamics-based Ising machine atp, such as simulated bifurcation and its variants, the approach acts as follows.
Firstly, the approach classifies the nodes into trapped and swing nodes. This is done by comparing |xp ,i | with the cutoff for all nodes i = 1, 2, ..., n. Then, the state xp is binarized by mapping swing nodes to zeros and trapped nodes xp ,j to σ j ∈ {+1, −1} based on the signs of the state values. This is motivated by the "freeze-out" of trapped nodes, in the sense that their magnitudes are sufficiently large such that further evolution of the system does not change the signs anymore.
Next, the approach targets to decide the remaining variables xp ,i , i ∈ V 0 , based on the trapped nodes. Swing nodes are randomly sequentially assigned binary values {+1, −1} based on the fields due to their neighbors, in accordance to the following two steps.
Firstly, every variable xp ,i , i ∈ V 0 is set to sign(y i ), with y i = n j=1,j =i G ij σ j by checking the stability of the configuration under the spin flip of node i (see Methods section). Secondly, the same stability criterion is applied to all nodes j = 1, 2, ..., n, leading to a stabilization of the whole configuration. Values of σ k are flipped to −σ k if the stability criterion σ k = sign(y i ) is violated. The updating order for both steps is randomized.

Supplementary Note 8 -Performance analysis of CIM with trapping-andcorrection (TAC) approach
In this section we demonstrate the performance of CIM with the proposed trapping-and-correction (TAC) approach on ten Max-cut instances, g05 100.0-g05 100.9 (unweighted graphs with n = 100 nodes with edge probability 0.5), from the Biq Mac Library [5] (with known solutions)complementary to the main text.
It can be observed that the mean Ising energy (over 100 trials) of the proposed approach always outperforms CIM (solid lines) and the top 5% results of the approach is clearly below both average and top 5% results of CIM.
Supplementary Figure 8: The figures compare average (solid lines) and top 5% quantile (dashed lines) Ising energies of ten Max-cut instances with n = 100 variables (100 random trials each). The grey lines are standard CIM, the blue lines are CIM with TAC approach. Ground state is marked as "min" and depicted by solid straight line.

Supplementary Note 9 -Performance analysis of simulated bifurcation and its variants with trapping-and-correction (TAC) approach
In this section, we demonstrate the performance of the TAC approach applied to other general Ising dynamics beyond CIM, such as simulated bifurcation (SB) and its variants, ballistic simulated bifurcation (bSB) and discrete simulated bifurcation (dSB). The size of test graphs, fully connected with random weights (±1), are increased to n = 2000 to examine possible performance variation in larger problem instances.
Supplementary Figure 11: Acceleration of discrete simulated bifurcation (dSB) through the trapping-andcorrection (TAC) approach. Here the Ising energies obtained by the SB algorithm is compared against those by TAC-dSB for random graphs of up to n = 2000 nodes (the graphs are fully connected with the weight of each edge drawn randomly from {−1, +1}). The solid curve and the dashed curve indicate the mean and top 5% quantile of the Ising energy respectively. The histograms illustrate the distribution of the obtained Ising energy at different steps.