A method to suppress local minima for symmetrical DOPO networks
AA method to suppress local minima for symmetrical DOPOnetworks
Seiya Amoh ∗ Daisuke Ito † Tetsushi Ueta ∗ Abstract
Coherent Ising machine (CIM) implemented by degenerate optical parametric oscillator (DOPO)networks can solve some combinatorial optimization problems. However, when the network structure hasa certain type of symmetry, optimal solutions are not always detected since the search process may betrapped by local minima. In addition, a uniform pump rate for DOPOs in the conventional operationcannot overcome this problem. In this paper proposes a method to avoid trapping of the local minima byapplying a control input in a pump rate of an appropriate node. This controller breaks the symmetricalproperty and causes to change the bifurcation structure temporarily, then it guides transient responsesinto the global minima. We show several numerical simulation results.
In recent years, the combinatorial optimization problems appears in various fields.There are a lot of algorithms to tackle them, but many Non-deterministic polynomial hard (NP-hard)problems require exponentially increasing computation time as the problem becomes larger.The Ising model is a mathematical formulation mimicking the ferromagnetism dynamics, and the searchingproblem of its ground state (namely the Ising problem[1]) can encode MAX-CUT problem[2, 3] for 3D networkgraphs. The solution obtained from the coherent Ising machine (CIM) takes a good value as an approximatesolution of the Ising problem. It has been demonstrated that Ising problems is solved by DOPO networkswhich is one of the implementations of CIM[4]. In the network, each DOPO corresponds to the Ising spinin the problem, and the network state is obtained as a solution when the pump rate exceeds the oscillationthreshold.Wang, et al. [5] provided theoretical background of the DOPO and showed the performance theoreticaland numerically, i.e., the coupled DOPOs can seek the optimal solution as a ground state of the Ising model.In most cases, the DOPO network gets the ground state by gradually pumping all oscillators. However, thesolution of the DOPO network is not always optimal[6]. In the case of gradually pumping DOPOs, localminima, which are not the ground state, appear earlier than the optimal solution, therefore they disturbfinding the ground state[7, 8]. As a result, the performance of the CIM decreases.Ito investigated bifurcation problems in a network configured by eight DOPOs[9]. By choosing a pumprate of a specific node and a unified pump rates of other nodes as a pair of parameters, several bifurcationdiagrams composed by super-critical and sub-critical pitchfork bifurcations and tangent bifurcations havebeen obtained. The analysis of their structures reveals that the variation of the unified pump rate tendsto be confined into a valley sectioned by bifurcation curves. Inside of the valley, the optimal solutions areconcealed by other local minima. This situation is caused by pitchfork bifurcations due to a symmetricalproperty, thus we guess that the performance of the CIM may be improved by breaking a symmetric propertyon purpose. In this paper, we propose a method to suppress local minima by adjusting not the unified pumprate but the specific one. This control perturbs the symmetrical property of the DOPO network temporarilyand guides the trajectory into the global solution. We show simulation results for several kinds of symmetricDOPO networks. ∗ Tokushima University, 2–1 Minami-Josanjima, Tokushima, 770-8501 Japan † Gifu University, 1–1 Yanado, Gifu, 501-1193 Japan a r X i v : . [ n li n . AO ] J u l DOPO networks
A DOPO system is described by the following c -number Langevin equations[4, 5]: dc j dt = ( p − c j − s j )) c j + N (cid:88) l =1 ,l (cid:54) = j ξ jl c l ,ds j dt = ( p − − ( c j − s j )) s j + N (cid:88) l =1 ,l (cid:54) = j ξ jl s l ,j = 1 , , . . . , N (1)where c and s are the normalized in-phase and quadrature-phase components, respectively. p is the pumprate, ξ is the coupling coefficient of DOPO network. A detailed description of the quantum models is shownin Refs.[10, 11]. In the numerical simulations of Eq. (1), we experimentally observe that s tend to be vanishedas time goes on. Thus the system can be simplified as follows: dc j dt = ( p − c j − c j + N (cid:88) l =1 ,l (cid:54) = j ξ jl c l . (2)If we define the state vector c and edge weight matrix Ξ as: c = c c ... c N , Ξ = ξ . . . ξ N ξ ξ N ... . . . ξ N ξ N , (3)where, ξ jl is a coupling coefficient between nodes j and l. If node i and j is coupled, ξ ij = ξ ji = ξ , but if not, ξ ij = ξ ji = 0. Then Eq. 1 can be derived as Eq. (4). d c dt = f ( c ) + Ξ c , (4)where f ( c ) = ( f ( c ) , f ( c ) , . . . , f ( c N )) T . In Eq. (2), the Jacobian matrix is derived as: J ( c ) = ∂ f ∂ c + 1 ξ Ξ . (5)MAX-CUT problem is a problem that divides nodes of a graph into two groups and maximizes the numberof edges cut. The Ising Hamiltonian can be applied to the number of cuts in the MAX-CUT problem[3]. Thenumber of cuts when describing the MAX-CUT problem in CIM using DOPO networks is defined as:cut(Ξ , χ ) = − N (cid:88) j =1 N (cid:88) l =1 ,l (cid:54) = j J jl χ j χ l , (6)where, χ j = c j / | c j | ., J jl is the generic element of the coupling matrix which defined as J jl = 0 if ξ jl = 0 ,otherwise J jl = 1. The relaxation function in Eq. (7) is the smallest value when the DOPO network is in theoptimal solution. η ( c ) = N (cid:88) j =1 ( p − c j −
1) = − N (cid:88) j =1 N (cid:88) l =1 ,l (cid:54) = j ξ jl c l c j . (7)Figure 1 gives MAX-CUT problems and their optimal solutions treated in this paper. The cut is de-termined by the sign of the solution c in Eq. (2), e.g., c , c , c , c take positive value and c , c , c , c takenegative value at Fig. (1)(a). Note that both graphs contain certain classes of symmetry and pitchforkbifurcations are the major bifurcation phenomena of the DOPO networks[12].2
23 456 78 (a) N = 8 (b) N = 10Figure 1: MAX-CUT problems and one of their optimal solutions. The number of edges which the dashedline intersects give the maximum for each graph. We may obtain candidates of solutions of the problem as equilibrium points of Eq. (1). From appropriateinitial values, the solution may fall into one of stable equilibrium points according to both global solutionsand local minima after a transition state. To compute equilibrium points, we apply Newton’s method byletting the left-hand of Eq. (2) be zero, with providing initial points by a brute-force strategy with Mersennetwista[13]. Figure 2(a) and (b) are stable equilibria of Eq. (2) corresponding to Fig. 1 (a) and (b), respectively.There are many unstable equilibria from the pitch-fork bifurcation around these attractors, but we do notvisualize them. Red and blue points are local minima and optimal solutions, respectively. The pump rate is p = 0 .
92 both of them. -1 -0.5 0 c c -1-0.510.50 -1 -0.5 0 c c -1-0.510.50 (a) (b)Figure 2: Stable equilibria. Blue points are the optimal solutions. (a): corresponding to Fig. 1(a), (b):corresponding to Fig. 1(b). Both are obtained with p = 0 . p from zero to one [14] comparedwith the fixed pump rate. That is, in the hardware implementation, we practically increase the pump ratefrom a low value, and gradually increase it near one monotonically.3 .1 Bifurcations of equilibria with pump rates Figure 3 shows the one-parameter bifurcation diagram for the case N = 8 corresponding to Fig. 1(a) andthis incremental pump rate process tracks this diagram from left to right. At any value of p , several randominitial values are scattered to find the equilibrium point.Gray points show the origin stable point corresponding to no cut. As p increases, the system receivesa pitch-fork bifurcation. Therefore, two symmetrical equilibrium points occur near the origin. In this case,they are corresponding to local minima. If p increases further, and at p = 0 . . < p <
1. There is no parameter set at which only the optimalsolution occurs. With this incremental pump rate process, local minima as sub-optimal solution is detectedat first, and the state is trapped to the same solution against further increment of the pump rate. It causesthe performance deterioration of the machine. − − − c p Figure 3: One-parameter bifurcation diagram ( N = 8).Let p j be an individual pump rate according to the node j and assume it is adjustable. Also assume thatand other pump rates are fixed. Figure 4 shows an incremental process of p for Fig. 1(a). The verticalaxis shows the amplitude of c . We increase p from 0.9 to 1 gradually. Blue and red point show optimalsolutions and local minima, respectively. The red solution meets a tangent bifurcation ( G ) near 0.99, andthen it disappears. Thus, there is parameter sets 0 . < p that gives only optimal solutions by adjustingthe pump rate individually.Now we investigate a relationship between optimal solutions and the pump rate. As mentioned above,not only global minima but also local minima are generated by pitchfork bifurcations induced from thesymmetrical configuration of the network. If the system is trapped into local minima, the value of therelaxation Eq. (7) becomes high since this stable situation is given by the forcing energy frustration to a partof nodes[6]. While, the precede work[9] suggested that a possibility of the control method to avoid trappingto local minima by adjusting pump rates of some DOPOs. This suppresses some bifurcation phenomena bybreaking the symmetry of the dynamical system on purpose. It can be called as a “bifurcation avoidancecontrol” so to speak. Figure 5 shows the p - p j bifurcation diagram at N = 8, where G and pf shows thetangent bifurcation and pitchfork bifurcation, respectively. We pick up a pump rate p j for one of DOPOs,and suppose a uniform pump rate p for other DOPOs. Only local minima occur in (i), both local minimaand the optimal solutions are mixed in (ii), and only optimal solutions occur in (iii).On the gray diagonal line in Fig. 5, the pump rate set does not reach (iii). However, controlling p independently, as indicated by the red arrow, the pump rate set reaches (iii).4 − − c p G Figure 4: c - p bifurcation diagram ( N = 8). p p pf (i) (iii) C on t r o lli ng (ii) pfG Figure 5: p - p j bifurcation diagram. In the Sect. 3.1, we recognize that there is the pump rate set in which only the optimal solution can beobtained by setting the pump rate p j independently. Although, with the method, it is necessary to specifythe amount to raise the pump rate in advance. The pump rate set which is only the optimal solution can beobtained from the relation with p j and c j . The frustration at a local minima imposed on nodes with small | c j | . To increase the | c j | , increase the pump rate p j of a node with the smallest | c j | . Therefore, consider afunction that automatically adjusts p j according to the size of | c j | , and replace p j with the pump rate ofthe node whose absolute value is small. Assume that a variable range of the pump rate is from zero to one,then e −| c j | can be adopted as a reasonable adjustment control value. As a result, the pump rate p j can beregarded as a function depending on the state variable, and the parameter can be set automatically.We propose a dynamic pump rate p j described as follows: p j = βe −| c j | , (8)where, β is gain of control. In this paper we set β = 1. When c j takes a smaller value, p j inflates. Figure 6is the change of p j . In the case of Fig. 5, p is controlled since | c | is the smallest. p rises according to the5xponential function of Fig. 6. As a result, the pump rate set moves to (iii). − − − − p j c j Figure 6: The example graph of Eq. (8), where β = 1.Figure 7 depicts the procedure of the control. When the calculation starts and the system becomes stablenear the equilibrium point, find the index j of the node with the smallest | c j | . Then, replace the pump rateof the j -th node with e −| c j | . The number of replacing pump rate is not specified in particular, excessivecontrol causes the performance degradation, so it is necessary to add control while monitoring the value ofthe relaxation function. In the following trials, we apply this procedure once. Replace p j with The solution is nearthe equilibrium pointFind j where is c j |Start NOYES
Relaxation function rises the smallest |End
NOYES e -| c j | Figure 7: The flowchart of the control.
Figure 8 shows the controlling result for N = 8 DOPO networks shown in Fig. 1(a) with the couplingcoefficient ξ = − .
1. We fix the unified pump rate as p = 0 .
92. As far as we checked, | c j | for all equilibrium6oints is less than unity, thus we put β = 1. The system runs without the controlling until t = 400, thenit converges the local minima, and has distribution of c − with large deviation. In the case of the trial inFig. 8, the node with j = 4 is the smallest | c j | , so p is selected as the control target. The controlling startsafter that, and the pump rate p is derived from Eq. (8). It is confirmed that the sign inversion of c occursafter t = 400, and the number of cuts is improved from 12 to 16 by this sign inversion. A red arrow in Fig. 5represents the change of p in the parameter plane. Note that, the local minima exist in the (i) region, theoptimal and local minima coexist in the (ii) region. In the (iii) region, the DOPO has only optimal solutions.The parameter set moves from the blue region to the red region by the proposed controller through thetangent bifurcation. As a result, our controller can avoid the local minima. By 100 trials with random initialvalues, our controller improve the probability of finding the optimal solution from 84% to 100%. -0.50 0.5 c - w/o controlling w/ controlling c sign inversion00.51 p
0 200 400 600 800 t
12 [cuts] 16 [cuts]
Figure 8: Control response for the N = 8 network. The controller adjusts the pump rate p .The control response for the N = 10 DOPO network is shown as Fig. 9. Here, the network structure isshown in Fig. 1(b), where the coupling coefficient ξ = − . p = 0 .
92. Wechoose the adjustable pump rate as p . Similar to Fig. 8, the control improves the number of cuts from 18to 26. -0.50 0.5 c - w/o controlling w/ controlling sign inversion p
0 200 400 600 800 t
18 [cuts] 26 [cuts]
Figure 9: Control response for the N = 10 network. The controller adjusts the pump rate p .Figure 10 shows the relationship between eight nodes on a steady-state with/without the proposed con-troller corresponding to the simulation results in Fig. 8. Note that, the size of each node is scaled by theabsolute values of c j of the node j . When the network converges to local minima solutions without controlling,the difference between maxima and minima circles becomes large; therefore node five and eight indirectlyfrustrate node one. Thus, they have | c j | values larger than other nodes. However, the difference in the size7f the nodes is relieved because the network converges to the global optimal solution from the local minimaby the proposed controller.Furthermore, Fig. 12 shows the transition of the relaxation function in each graph. The change of thenode size in the graph of Fig. 10 is evaluated by Fig. 12(a). The η ( c ) of Fig. 10 decreases at t = 400 due tothe addition of control. Fig. 12(c) and (d) are N = 16 and N = 20 nodes graph which are Fig. 11. Similarto Fig. 1(a) and (b), Fig. 11 also have symmetry. Initial pump rate of N = 16 and N = 20 are p = 0 . N = 16 and N = 20 improves from 28 to 32 and 36 to 40, respectively. (a) (b)Figure 10: (a): before controlling, (b): after controlling. (a) N = 16 (b) N = 20Figure 11: The graph of N = 16 and N = 20.Table 1 summarizes the probability of the optimal state with and without the proposed controller. Forthe four networks, the proposed controller can suppress the local minima, and can improve the probabilityof the optimal solutions.Table 1: Probability of the optimal solution for w/o and w/ proposed controller. The number of trials is 100.The initial state is given by random values.number of nodes w/o control w/ control8 84% 100%10 78% 100%16 92% 97%20 88% 100%8 − − − − η ( c ) t
12 [cuts] 16 [cuts] − − − − − η ( c ) t
28 [cuts] 32 [cuts] (a) (b) − − − − − η ( c ) t
28 [cuts] 32 [cuts] − − − − − η ( c ) t
36 [cuts] 40 [cuts] (c) (d)Figure 12: Transitions of the relaxation function.(a): N = 8, (b): N = 10, (c): N = 16, (d): N = 20. We propose a method to avoid local minima for symmetrical DOPO networks. We design new controllingscheme for the improvement of the DOPO network as the MAX-CUT solver. We summarize the proposalmethod and the results to following items.1. Conventionally, the pump rate was set uniformly, although we recognize that the local minima disap-peared due to tangent bifurcation when the pump rate was set individually at each node.2. There is a pump rate set that can obtain only the optimal solution in the p j - p bifurcation diagram ofthe DOPO network.3. From the bifurcation diagram and the structure in the DOPO network, we proposed a controllingmethod that automatically sets the pump rate by the exponential function.4. Our controller improves the probability of success for optimal solution search in the four symmetricalDOPO networks.Next, the future tasks are follows:.1. In this paper, the relaxation function was monitored for the number of replacing pump rate, althoughwe need to find out how to determine the most effective number of replacing.2. We try to systematic explanations of the performance of proposed controller based on the bifurcationtheory, however since there is a non-differentiable absolute value function in the control term, thebifurcation analysis becomes difficult. This problem can be avoided by replacing e −| c j | with a quadraticfunction such as 1 − c j , but the performance comparison between two control terms is required.This paper is under article submission and will be published on Nonlinear Theory and Its Applications,IEICE(Vol.E11-N, No.4, Oct. 2020). All authors agree to post this preprint to arXiv.9 eferences [1] E. Ising, “Beitrag zur theorie des ferromagnetismus,” Zeitschrift f¨ur Physik , vol. 31, no. 1, pp. 253–258,1925.[2] R. M. Karp, “Reducibility among combinatorial problems,” in
Complexity of Computer Computations ,pp. 85–103, Boston, MA: Springer US, 1972.[3] A. Galluccio, M. Loebl, and J. Vondr´ak, “New algorithm for the Ising problem: Partition function forfinite lattice graphs,”
Physical Review Letter , vol. 84, no. 26, pp. 5924–5927, 2000.[4] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, “Coherent Ising machine based ondegenerate optical parametric oscillators,”
Physical Review A , vol. 88, no. 6, p. 063853 (9 pages), 2013.[5] Z. Wang,
Coherent computation in Degenerate Optical Parametric Oscillators . PhD thesis, StanfordUniversity, Stanford University, Stanford, CA 94305 USA, 2015.[6] Y. Haribara, S. Utsunomiya, and Y. Yamamoto, “A coherent Ising machine for MAX-CUT problems:Performance evaluation against semidefinite programming and simulated annealing,” in
Principles andMethods of Quantum Information Technologies , vol. 911 of
Lecture Notes in Physics , ch. 12, pp. 251–262,Springer, Japan, 2016.[7] K. Takata, S. Utsunomiya, and Y. Yamamoto, “Transient time of an Ising machine based on injection-locked laser network,”
New Journal of Physics , vol. 14, p. 013052, 2012.[8] K. Takata, A. Marandi, R. Hamerly, Y. Haribara, D. Maruo, S. Tamate, H. Sakaguchi, S. Utsunomiya,and Y. Yamamoto, “A 16-bit coherent Ising machine for one-dimensional ring and cubic graph problems,”
Scientific Reports , vol. 6, p. 34089, 2016.[9] D. Ito, T. Ueta, and K. Aihara, “Bifurcation analysis of eight coupled degenerate optical parametricoscillators,”
Physica D: Nonlinear Phenomena , vol. 372, pp. 22–30, 2018.[10] T. Shoji, K. Aihara, and Y. Yamamoto, “Quantum model for coherent Ising machines: Stochasticdifferential equations with replicator dynamics,”
Phys. Rev. A , vol. 96, no. 5, p. 053833, 2017.[11] A. Yamamura, K. Aihara, and Y. Yamamoto, “Quantum model for coherent Ising machines: Discrete-time measurement feedback formulation,”
Phys. Rev. A , vol. 96, no. 5, p. 053834, 2017.[12] T. Ueta, H. Miyazaki, T. Kousaka, and H. Kawakami, “Bifurcation and chaos in coupled BVP oscilla-tors,”
International Journal of Bifurcation and Chaos , vol. 14, no. 04, pp. 1305–1324, 2004.[13] M. Matsumoto and T. Nishimura, “Mersenne twister: A 623-dimensionally equidistributed uniformpseudo-random number generator,”
ACM Transactions on Modeling and Computer Simulation , vol. 8,no. 1, pp. 3–30, 1998.[14] K. Takata and Y. Yamamoto, “Data search by a coherent Ising machine based on an injection-lockedlaser network with gradual pumping or coupling,”