Discrete Polynomial Optimization with Coherent Networks of Condensates and Complex Coupling Switching
DDiscrete Polynomial Optimization with Coherent Networks of Condensates andComplex Coupling Switching
Nikita Stroev and Natalia G. Berloff , ∗ Skolkovo Institute of Science and Technology, Bolshoy Boulevard 30, bld.1, Moscow, 121205, Russian Federation and Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Cambridge CB3 0WA, United Kingdom
Gain-dissipative platforms consisting of lasers, optical parametric oscillators and nonequilibriumcondensates operating at the condensation/coherence threshold have been recently proposed as effi-cient analog simulators of 2-local spin Hamiltonians with continuous or discrete degrees of freedom.We show that nonequilibrium condensates above the threshold arranged in an interacting networkmay realise k -local Hamiltonians with k > Recently much effort has been devoted to the devel-opment of various technological platforms that act asquantum or classical analog simulators aimed at solvingcertain classes of hard classical optimization problems[1–7]. It is expected that these kinds of platforms wouldhelp to efficiently solve many tasks of significant compu-tational complexity, ranging from modelling microscopiceffects and processes like the behavior of electrons in com-plex materials [8, 9] and finding the ground state of spinglasses [10], to the applied combinatorial optimizationproblems [11]. Large scale computational problems ofthis type are hard for classical von Neumann architec-ture which suggests looking for fully analog or hybriddigital/analog/quantum devices that can find a solutionfaster or find a better solution in a fixed time.Nonequilibrium condensates, optical parametric oscil-lators, lasers, and memristors have been considered asannealing-inspired accelerators and demonstrated suc-cesses in finding the ground state of spin Hamiltonianswith continuous or discrete variables [2, 6, 12, 13]. Inparticular, the Coherent Ising Machine has been shownto significantly outperform classical simulated annealingin terms of both accuracy and computation time to effi-ciently solve Max-CUT problems [2] and has shown bet-ter scalability than the quantum annealers [14]. Mem-ristors with massively parallel operations performed ina dense crossbar array were shown to be able to solveNP-hard Max-CUT problems predicting over four ordersof magnitude advantage over digital and optical hard-ware [13]. Integrated photonic circuits that use self-phasemodulation in two microring resonators were shown toact as an optical coherent Ising machine [15, 16]. Thelattices of exciton-polariton condensates were shown toefficiently simulate the XY Hamiltonian when operatingat the condensation threshold [12, 17] In all these sys- ∗ correspondence address: N.G.Berloff@damtp.cam.ac.uk tems, discrete Ising ’spins’ or continuous XY ’spins’ areencoded in individual phase modes of the nonlinear net-works. An optimization problem of interest is mappedinto the connectivity matrix of the spin network with thetask of finding its ground state, which can be related tofinding the ’maximum occupancy’ of the collective su-permode of the underlying network, as a system specificgain mechanism is continuously increased to reach thecoherence threshold [12, 18].The development of technological platforms thatpromise to offer a significant time or power consump-tion improvements in solving hard optimization problemsgoes hand in hand with annealing inspired optimization,when the physical principle of the device is used to for-mulate the classical algorithms [19]. The recent exam-ples of which include Simulated Bifurcation algorithmsinspired by quantum adiabatic optimization using a non-linear oscillator network [20], destabilization of local min-ima based on degenerate parametric oscillator networks[21], parallel tempering Monte Carlo [22], and the gain-dissipative algorithm based on operation of the polaritongraph simulator at the condensation threshold [23].The focus of all these technological and inspired imple-mentations of the annealer-based optimization has beenon QUBO, however, there is a large class of optimiza-tion problems — the higher order polynomial binary op-timization (HOBO) – that are more naturally encodedby the k -local Hamiltonians [24, 25]. HOBO is concernedwith optimizing a (high degree) multivariate polynomialfunction in binary variables. Our basic model is to maxi-mize or minimize a k -th degree polynomial function f ( x )where x = ( x , ..., x i , ..., x N ), x i ∈ {− , +1 } . The ex-amples of HOBO are ubiquitous from Hypergraph max-covering problem to Frobenius and ”market split” prob-lems [25]. HOBO is a fundamental problem in integerprogramming and is also known as Fourier support graphproblem. Any HOBO can be mapped into the QUBO[26], however, the overhead in the number of nodes be-comes prohibitive in an actual technological platform, so a r X i v : . [ c ond - m a t . d i s - nn ] M a y it is important to consider ways to solve HOBO directly.The purpose of this article is three-fold. First, we showthat Ising machines based on nonequilibrium condensatescan be used to address 4-local HOBO when operating above the threshold. Secondly, inspired by the operationof the networks of nonequilibrium condensates we pro-pose a new optimization algorithm for solving HOBO ofarbitrary degree. Finally, we show that another physics-inspired method of turning on and off the complex cou-pling between the nonlinear condensates greatly enhancesthe search for global minimum. Polynomial optimization with coherent networks.
Theoptimization problem studied in this paper ismin x ∈{− , +1 } N − (cid:88) Ω A ki ,i , ··· ,i k x i x i · · · x i k , (1)where Ω = { i j : 1 ≤ i ≤ i ≤ ... ≤ i k ≤ N } and A k isthe super-symmetric tensor of degree k . To formulate thegain-dissipative platform that reaches the ground state ofHOBO by finding the ’maximum occupancy’ collectivesupermode of the underlying network of nonequilibriumcondensates we consider the mean-field equations thatgovern such a network based on the Ginzburg-Landauequation [27, 28]. This is a universal driven-dissipativeequation that describes the behaviour of systems in thevicinity of a symmetry–breaking instability and has beenused to describe lasers, thermal convection, nematic liq-uid crystals, and various non-equilibrium condensates[29, 30]. When derived asymptotically from a genericlaser model given by Maxwell-Bloch equations it has asaturable nonlinearity and can be written asi ∂ψ∂t = −∇ ψ + ˜ U | ψ | ψ + i (cid:18) P ( r , t )1 + b | ψ | − γ c (cid:19) ψ, (2)where ψ ( r , t ) is the wavefunction of the system, ˜ U is thestrength of the delta-function interaction potential, γ c is the rate of linear losses, b parametrizes the effectivestrength of nonlinear losses, P ( r , t ) describes the gainmechanism that adds particles to the system. It was ex-perimentally demonstrated [12] that when pumped at thecondensation threshold, freely expanding optically im-printed polariton condensates arranged in a lattice mayachieve a steady state with condensate phases realisingthe minimum of the XY Hamiltonian. In this frame-work, the coupling strengths between condensates de-pend of the system parameters, pumping intensity andshape and on the lattice geometry [31]. We shall assumethat P ( r , t ) adds particles in N spatial locations cen-tered at r i , i = 1 , ..., N , so that P ( r , t ) = (cid:80) i f i ( t ) p i ( r ),where f i is the time-dependent part of the pumping at r = r i and p i ( r ) ≡ p ( | r − r i | ) is a given spatially lo-calised pumping profile, that creates the condensate witha wavefunction φ i ( r ) ≡ φ ( | r − r i | ) centred at r = r i and normalized so that 2 π (cid:82) | φ ( r ) | r dr = 1. In writingEq. (2) we let (cid:126) = 1 and m = 1 /
2. If the distances be-tween the neighbouring condensates are larger than the width of p ( r ), we employ the tight binding approxima-tion and write the wavefunction of the system as a linearsuperposition of the wavefunctions of individual coherentcenters ψ ( r , t ) ≈ (cid:80) Ni =1 a i ( t ) φ i ( r ), where a i ( t ) is the time-dependent complex amplitude networks [23, 32]. We ex-pand the first term in the brackets of Eq. (2) in Taylorseries, substitute the expressions for P and ψ , multiplyby φ ∗ j for j = 1 , ..., N and eliminate the spatial degreesof freedom by integrating in the entire space to obtain N equations of the form d Ψ i dt = Ψ i (cid:16) γ i − ( σ i + i U ) | Ψ i | (cid:17) + (cid:88) j,j (cid:54) = i J ij Ψ j + (cid:88) (cid:104) j,k,l (cid:105) Q ijkl Ψ j Ψ k Ψ ∗ l . (3)In writing Eq. (3) we used the following notations: γ i = f i (cid:82) p | φ | d r − γ c , σ i = f i (cid:82) p | φ | d r , U = ˜ U (cid:82) | φ | d r , Q ijkl = − b (cid:80) m ∈{ i,j,k,l } f m (cid:82) p m φ k φ ∗ l φ i φ ∗ j d r , and (cid:104) i, j, k (cid:105) denotes the combinations of { j, k, l } that exclude j = k = l = i . To the leading order we kept theterms J ij = Re [ f i (cid:82) p i φ i φ ∗ j d r + f j (cid:82) p i φ j φ ∗ i d r ] for j (cid:54) = i, and denoted Ψ i = a i exp(i tRe (cid:82) φ ∗ ∇ φ d r ). Werewrite Eq. (3) in terms of the number densities ρ i andphases θ i , Ψ i = √ ρ i exp[i θ i ] , ˙ ρ i ( t ) = ( γ i − σ i ρ i ) ρ i + (cid:80) j,j (cid:54) = i J ij √ ρ i ρ j cos θ ij + (cid:80) (cid:104) j,k,l (cid:105) Q ijkl √ ρ i ρ j ρ k ρ l cos θ ijkl , and ˙ θ i ( t ) = − U ρ i − (cid:80) j,j (cid:54) = i J ij (cid:112) ρ j /ρ i sin θ ij − (cid:80) (cid:104) j,k,l (cid:105) Q ijkl (cid:112) ρ j ρ k ρ l /ρ i sin θ ijkl , where θ ij = θ i − θ j and θ ijkl = θ i + θ l − θ k − θ j .The higher order terms affect the states even in thesimplest configuration of two identical oscillators pumpedwith γ i = γ for which the occupancy ρ = ρ = ρ at thefixed point reads ρ = [ γ + J cos ∆ θ + ˜ Q cos(2∆ θ )] /σ where ∆ θ = θ − θ and ˜ Q = ρ Q. By choosingthe minimum pumping γ to reach the required occu-pancy, we minimize the Hamiltonian H two = − J cos ∆ θ − ˜ Q cos(2∆ θ ) while the equation on θ i describes the gradi-ent descent to the local minimum of H two with the dy-namics of the higher order Kuramoto oscillators. If ˜ Q isnegligible (close to condensation threshold), we have theminimization of the XY Hamiltonian, so ∆ θ = 0 or π if J >
J < Q is present but has the same sign as J . How-ever, a different phase difference is realised when J ˜ Q < | ˜ Q/J | ≥ , namely ∆ θ = arccos( − J/ Q ).In the example of two oscillators the stationary statewith equal occupancy of the nodes is always reached.However, in a more general system with many oscilla-tors, unless the oscillatory network is highly symmetric(all oscillators have equal number of connections of thesame strength with other oscillators) the systems breaksinto subsystems characterised by different frequencies.To guarantee the full synchronisation of the network weneed to choose the injection rates in such a way thatall oscillators have the same occupancy ρ th [17]. For in-stance, this can be achieved by adjusting the pumpingrates dynamically, depending on the occupancy of the i − th oscillator at time t :˙ γ i = (cid:15) ( ρ th − ρ i ) , (4)where the parameter (cid:15) characterizes how fast γ i adjusts tochanges in ρ i . Aiming at the algorithmic implementation,we will focus on tensors of the same order k , and leave theproblems with mixed order tensors for future work. Incase of polariton condensates, this means that the fourthorder tensors dominate the dynamics of the second-order terms. With appropriate density adjustments, asdescribed above, the system of N oscillators will al-ways synchronize and achieve the stationary minimumof the Hamiltonian H = − (cid:80) i (cid:80) (cid:104) j,k,l (cid:105) ˜ Q ijkl cos θ ijkl withsuper-symmetric tensor of k = 4.To replace the minimization in the space of continuousspins with binary states, one could combine nonresonantpumping with resonant at twice the frequency of the con-densate which introduces the terms proportional to ψ ∗ (Ψ ∗ i ) to the right-hand side of Eq. (2) (Eq. (3)) similarlyto k = 2 case [17]. Resonant and nonresonant excita-tions have been previously combined in experiments onpolariton condensates using chemical etching across thesample allowing resonant excitation from the back sideof the cavity [33]. Such combination of resonant andnonresonant excitations would lead to the realisation ofthe minimum of the k -local Hamiltonians, so to solvingEq. (1) with the binary spins x i = cos θ i , where θ i arelimited to 0 and π . However, in contrast with k = 2 case,the phase projections on the binary states for k > i j and Ψ ∗ i k present inthe tensor form as shown below, so the presence of theadditional resonant field in not necessary. Physics-Inspired Optimization.
The principle of coher-ence formation at a minimum of a spin Hamiltonian for-mulated above inspires an efficient algorithm for findingthe global minimum of HOBO. For this, we extend andsimplify Eq. (3) to capture the mechanism of relaxationto the minimum of the HOBO but without the necessityto capture full physics of the actual system. The mini-mum of HOBO for N binary variables can be found bynumerical integration of 2 N equations d Ψ l dt = Ψ l ( γ l ( t ) − | Ψ l | ) + (cid:88) ¯Ω A ki ,i , ··· ,l, ··· i k Ψ i Ψ i · · · Ψ ∗ i k , (5)together with Eq. (4) where ¯Ω = Ω /l and theinitial values for pumping strength γ l ( t = 0) = − max ≤ l ≤ N (cid:80) ¯Ω | A ki ,i , ··· ,l, ··· i k | . At the fixed point, theimaginary part of Eq. (5) gives the set of linear equa-tions such that the l -th equation involves superposi-tion of sin( (cid:80) i j (cid:54) = { l,i k } θ i j − θ i k − θ l ) that has to beequal to zero. In general, the only way for the sys-tem to satisfy these equations is to bring all phases θ l to take on 0 or π . The total occupancy of the sys-tem at the fixed point is found from the real part of Eq. (4) and is equal to N ρ th , so that N ρ th = (cid:80) Nl =1 γ l + √ ρ th k − (cid:80) Ω A ki ,i , ··· ,i k cos( θ i ) cos( θ i ) · · · cos( θ i k ) . Ifwe set the process of raising the pumping from belowthat guarantees that (cid:80) Nl =1 γ l is the smallest possibleinjected intensity, then at the fixed point the systemfinds the global minimum of the k − local Hamiltonian H = − (cid:80) Ω A ki ,i , ··· ,i k cos( θ i ) cos( θ i ) · · · cos( θ i k ) , and,therefore, solves Eq. (1). We will refer to the Eqs. (5-4)as the Tensor Gain-Dissipative (TGD) method.To illustrate the behavior of the system we first con-sider a toy problem: the following 3-local Hamiltonian H test ( x ) = − x x x − x x x − x x x − x x x , with variables x i ∈ {± } , while Eq. (5) becomes ˙Ψ l =Ψ l ( γ l − | Ψ l | ) + (cid:80) (cid:104) j,k (cid:105) K ljk Ψ j Ψ ∗ k , and K is a tensor withnonzero entries 1 , , , . The Hamiltonian H test has 2 stationary points, among which there are three local min-ima: H = − , H = − , H = −
13 and the globalminimum H = −
15, that all can be accessed duringthe time evolution of the system. To understand thebasins of attraction for these stationary points we numer-ically integrate Eqs. (5-4) starting with initial conditionsΨ i ( t = 0) = exp[i θ i ] where the phases θ i ∈ [0 , π )are uniformly distributed [34]. Figure 1(a) depicts thestatistics of distribution of the stationary points reachedand indicates that the basins of local minima combinedare larger than that of the global minimum. To facilitatethe search for the global minimum the algorithm needs toallow for a possibility to explore the hyperspace until thelowest lying energy state is found. This can be achievedby adding a Langevin noise, which represents intrinsicvacuum fluctuations and classical noise, that shifts thetrajectory from its deterministic path while allowing itto stay below any local minima. Figure 1(b) depicts thestatistics of reaching local and global minima found bynumerical integration of Eqs. (5-4) using the same initialconditions as in Fig. 1(a) but with the white noise added.Decreasing the (cid:15) parameter allows to improve the possi-bility of reaching the global minima. Figre 1(d) showsone such trajectory as it approaches the global minimumof H test from below.With the growth in the number of variables and con-comitant growth of the system hyperspace local noisyperturbation of the trajectory may not be sufficient toreach the global minimum basin of attraction. Moti-vated by recent studies of heteroclinic networks [35] weintroduced heteroclinic orbits into our model by engi-neering time-dependent complex couplings into the net-work Eqs. (5-4). Complex couplings naturally appearin polariton model if the energy shift due to a noncon-densed reservoir R ( r , t ) is present in the system [36]. Theimaginary part of the coupling may destabilise the stablefixed point so that the system trajectory quickly leavesits neighborhood along the fastest direction. Includingthis switching dynamics into the system facilitates thesearch for the true global minimum by allowing fullerexploration of the phase space. Complex coupling switching . To implement the com-
FIG. 1: Success rates for achieving local and global minimaof H test using numerical simulations of Eqs. (5-4) using fullydeterministic integration without noise (a); with white noiseadded to the right hand side of Eq. (5) (b); using complexcoupling control as described in the main text (c). The in-dexes marked on the horizontal axis label local minima as1 , , H , H and H respectively. Theglobal minimum ( H ) is labeled 4. The insert (d) shows thetwo-dimensional projection ( θ = θ = 0) of the energy land-scape for H test with x i = cos θ i . The time evolution of thescaled total injected intensity ( (cid:80) γ i − Nρ th ) / plex coupling switching method (TGD+CC) on H test weturn two of the real coupling coefficient into the com-plex ones with a significant imaginary part as soon asthe system reaches a steady state. The system trajec-tory leaves the basin of attraction of that state and trav-els to a different part of the system hypercube [0 , π ] N .When the imaginary part of the coupling is turned offanother steady state will be found. By varying thecoupling elements to be switched, the time durationof the switching, and the amplitude of the imaginarypart while keeping the injected intensity low we allowthe system to efficiently search for a low energy mini-mum. In our test example, implementing the switch-ing of a coupling coefficients K and K accord-ing to K ( t ) = 8(1 + 4 i ) , K ( t ) = 4(1 − i ) , t ∈ [ t , t + 160] ∪ [ t , t + 160] ∪ [ t , t + 280] and keeping K ( t ) = 8 , K ( t ) = 4 otherwise allows every trajec-tory irrespective of its initial state to reach the globalminimum, see Fig. 1(c) . Here t , t , t are times at whichthe system settles to a steady state after switching off theimaginary part of the couplings [37]. Complex coupling switching for large N . We tested thecomplex couplings switching approach using the largescale simulations on 20 dense and 20 sparse random ten-sor sets K ijk of 3 d rank with 10 elements over 500 dif-ferent realisations. To implement TGD+CC method onlarge N, as soon as the system reaches the steady state FIG. 2: Success probability density distribution for 500 real-isations on the set of 20 dense (main) and sparse (inset) 3 d rank tensors K ijk with 10 elements for four different meth-ods discussed in the main text. To generate sparse tensors9/10th of all elements in dense tensors were set to zero. we randomly choose N/
50 of the coupling strengths K ijk and modify them by adding 3 iK ijk . As the system tra-jectory leaves the basin of attraction of the previous fixedpoint, we let return to original couplings and allow thesystem to relax to a new steady state. Keeping the totalinjected rate (cid:80) γ l small forces the trajectories to explorethe low energy states of the Hamiltonian until the trueglobal minimum is found [38].We compare the behaviour of the TGD and TGD+CCwith two popular network-based methods and show thatTGD+CC outperforms all of these methods. The firstmethod represents the network of analogue bistable units(NBU) in the presence of a double-well potential deriva-tive that forces the network elements x l to take on ± dx l dt = − hx l | x l | k − ( x l −
1) + (cid:88) ¯Ω A ki ,i , ··· ,l, ··· ,i k x i x i · · · x i k , (6)where x l ( t = 0) are randomly distributed real numbers,and h is a control parameter [39]. In comparison withthe usual k = 2 case [21], we balanced the degrees ofpolynomial between two term on the right-hand side ofEq. (6) by introducing | x l | k − factor. Another efficientsolver of Eq. (1) is given by a higher order Hopfield neuralnetworks [40]: dx l dt = − x l + (cid:88) ¯Ω A ki , ··· ,l, ··· ,i k tanh (cid:18) x i β (cid:19) · · · tanh (cid:18) x i k β (cid:19) , (7)where x l are real continuous variables and β is the scalingparameter [41].Figure 2 shows the results of large-scale numerical sim-ulations using different methods. TGD+CC consistentlyhas a better success probability of finding the global min-imum. Conclusions.
In this Letter we showed that latticesof nonequilibrium condensates when pumped above thethreshold may realise k -local Hamiltonians with non-trivial spin structures. We formulated system-inspiredmethod of computing the optimal solution of a largerange of HOBO problems. Finally, we introduced theconcept of computation via the mechanism of complexcoupling switching. Its combination with the tensor vari-ation of the gain-dissipative algorithm leads to an effi-cient way of finding the low energy states of the Hamil-tonian due to: (i) individual node gain control that allowsto explore low energy states of the Hamiltonian and guar- antees the achievement of the minimum, (ii) evolution inreal number space that allows to tunnel through func-tional barriers in discrete variables, (iii) Kuramoto net-works graduate decent close to threshold, and finally (iv)complex coupling switching that allows the trajectory toescape local minima. This approach offers a highly flex-ible new kind of computation that is inspired by andcompatible with many physical network realisations. Weenvision it becoming a part of a hybrid platform wherethe states of the network are fed to the physical devicefor the optimal performance.The authors acknowledge support from Huawei. [1] S. Utsunomiya, K. Takata, and Y. Yamamoto, “Mappingof ising models onto injection-locked laser systems,” Op-tics express , vol. 19, no. 19, pp. 18091–18108, 2011.[2] A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Ya-mamoto, “Network of time-multiplexed optical paramet-ric oscillators as a coherent ising machine,”
Nature Pho-tonics , vol. 8, no. 12, p. 937, 2014.[3] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Ya-mamoto, and H. Takesue, “Large-scale ising spin networkbased on degenerate optical parametric oscillators,”
Na-ture Photonics , vol. 10, no. 6, p. 415, 2016.[4] P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly,C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Ut-sunomiya, K. Aihara, et al. , “A fully programmable 100-spin coherent ising machine with all-to-all connections,”
Science , vol. 354, no. 6312, pp. 614–617, 2016.[5] Y. Takeda, S. Tamate, Y. Yamamoto, H. Takesue, T. In-agaki, and S. Utsunomiya, “Boltzmann sampling for anxy model using a non-degenerate optical parametric oscil-lator network,”
Quantum Science and Technology , vol. 3,no. 1, p. 014004, 2017.[6] M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson,“Observing geometric frustration with thousands of cou-pled lasers,”
Physical review letters , vol. 110, no. 18,p. 184102, 2013.[7] D. Dung, C. Kurtscheid, T. Damm, J. Schmitt,F. Vewinger, M. Weitz, and J. Klaers, “Variable poten-tials for thermalized light and coupled condensates,”
Na-ture Photonics , vol. 11, no. 9, p. 565, 2017.[8] I. Buluta, S. Ashhab, and F. Nori, “Natural and artificialatoms for quantum computation,”
Reports on Progress inPhysics , vol. 74, no. 10, p. 104401, 2011.[9] I. M. Georgescu, S. Ashhab, and F. Nori, “Quantumsimulation,”
Reviews of Modern Physics , vol. 86, no. 1,p. 153, 2014.[10] F. Barahona, “On the computational complexity of isingspin glass models,”
Journal of Physics A: Mathematicaland General , vol. 15, no. 10, p. 3241, 1982.[11] E. L. Lawler, J. K. Lenstra, A. R. Kan, D. B. Shmoys, et al. , The traveling salesman problem: a guided tourof combinatorial optimization , vol. 3. Wiley New York,1985.[12] N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D.T¨opfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis,“Realizing the classical xy hamiltonian in polariton sim- ulators,”
Nature materials , vol. 16, no. 11, p. 1120, 2017.[13] F. Cai, S. Kumar, T. Van Vaerenbergh, R. Liu, C. Li,S. Yu, Q. Xia, J. J. Yang, R. Beausoleil, W. Lu, et al. ,“Harnessing intrinsic noise in memristor hopfield neuralnetworks for combinatorial optimization,” arXiv preprintarXiv:1903.11194 , 2019.[14] R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli,A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba,T. Honjo, et al. , “Experimental investigation of perfor-mance differences between coherent ising machines anda quantum annealer,” arXiv preprint arXiv:1805.05217 ,2018.[15] D. Kielpinski, R. Bose, J. Pelc, T. Van Vaerenbergh,G. Mendoza, N. Tezak, and R. G. Beausoleil, “Infor-mation processing with large-scale optical integrated cir-cuits,” in , pp. 1–4, IEEE, 2016.[16] N. Tezak, T. Van Vaerenbergh, J. S. Pelc, G. J. Men-doza, D. Kielpinski, H. Mabuchi, and R. G. Beausoleil,“Integrated coherent ising machines based on self-phasemodulation in microring resonators,”
IEEE Journal ofSelected Topics in Quantum Electronics , vol. 26, no. 1,pp. 1–15, 2019.[17] K. P. Kalinin and N. G. Berloff, “Networks of non-equilibrium condensates for global optimization,”
NewJournal of Physics , vol. 20, no. 11, p. 113023, 2018.[18] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Ya-mamoto, “Coherent ising machine based on degener-ate optical parametric oscillators,”
Physical Review A ,vol. 88, no. 6, p. 063853, 2013.[19] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa,H. Tamura, and H. Katzgrabeer, “Physics-inspired op-timization for quadratic unconstrained problems using adigital annealer,”
Frontiers in Physics , vol. 7, p. 48, 2019.[20] H. Goto, K. Tatsumura, and A. R. Dixon, “Combinato-rial optimization by simulating adiabatic bifurcations innonlinear hamiltonian systems,”
Science advances , vol. 5,no. 4, p. eaav2372, 2019.[21] T. Leleu, Y. Yamamoto, P. L. McMahon, and K. Aihara,“Destabilization of local minima in analog spin systemsby correction of amplitude heterogeneity,”
Physical re-view letters , vol. 122, no. 4, p. 040607, 2019.[22] I. Rozada, M. Aramon, J. Machta, and H. G. Katz-graber, “Effects of setting the temperatures in the par-allel tempering monte carlo algorithm,” arXiv preprint arXiv:1907.03906 , 2019.[23] K. P. Kalinin and N. G. Berloff, “Simulating ising andpotts models and external fields with non-equilibriumcondensates,” arXiv preprint arXiv:1806.11457 , 2018.[24] S. He, Z. Li, and S. Zhang, “Approximation algorithmsfor discrete polynomial optimization,”
Journal of the Op-erations Research Society of China , vol. 1, no. 1, pp. 3–36,2013.[25] B. Jiang, Z. Li, and S. Zhang, “Approximation methodsfor complex polynomial optimization,”
ComputationalOptimization and Applications , vol. 59, no. 1-2, pp. 219–248, 2014.[26] G. De las Cuevas and T. S. Cubitt, “Simple univer-sal models capture all classical spin physics,”
Science ,vol. 351, no. 6278, pp. 1180–1183, 2016.[27] M. Wouters and I. Carusotto, “Excitations in a nonequi-librium bose-einstein condensate of exciton polaritons,”
Physical review letters , vol. 99, no. 14, p. 140402, 2007.[28] J. Keeling and N. G. Berloff, “Spontaneous rotating vor-tex lattices in a pumped decaying condensate,”
Physicalreview letters , vol. 100, no. 25, p. 250401, 2008.[29] I. Carusotto and C. Ciuti, “Quantum fluids of light,”
Reviews of Modern Physics , vol. 85, no. 1, p. 299, 2013.[30] J. Keeling and N. G. Berloff, “Exciton–polariton conden-sation,”
Contemporary Physics , vol. 52, no. 2, pp. 131–151, 2011.[31] P. G. Lagoudakis and N. G. Berloff, “A polariton graphsimulator,”
New Journal of Physics , vol. 19, no. 12,p. 125008, 2017.[32] A. Smerzi and A. Trombettoni, “Nonlinear tight-bindingapproximation for bose-einstein condensates in a lattice,”
Physical Review A , vol. 68, no. 2, p. 023613, 2003.[33] H. Ohadi, R. Gregory, T. Freegarde, Y. Rubo, A. Ka-vokin, N. G. Berloff, and P. Lagoudakis, “Nontrivialphase coupling in polariton multiplets,”
Physical ReviewX , vol. 6, no. 3, p. 031032, 2016.[34] Eqs. (5-4) were Euler integrated with (cid:15) = 0 . M , γ l ( t =0) = − . M , ρ th = 1 , with D = max(1 − ρ/ρ th ,
0) , dt =0 . T stop = 300, M = max ≤ l ≤ N (cid:80) ¯Ω | A kl,i ,i , ··· ,i k | .Tensor elements were scaled by 0 .
05 to slow the dynamicsdown.[35] F. S. Neves and M. Timme, “Computation by switchingin complex networks of states,”
Physical review letters ,vol. 109, no. 1, p. 018701, 2012.[36] K. P. Kalinin and N. G. Berloff, “Polaritonic network asa paradigm for dynamics of coupled oscillators,” arXivpreprint arXiv:1902.09142 , 2019.[37] The switching times are defined by analysing the transi-tions between the steady states labeled by 1 , , , K ( t ) = 8(1 + 4 i ) and K ( t ) = 4(1 − i )for t = T after the steady states are reached leads tothe following permutation of the states: a)(1 , , , → (2 , , ,
4) if T = 80; b)(1 , , , → (1 , , ,
2) if T = 160;c)(1 , , , → (4 , , ,
4) if T = 280 . Clearly the switch-ing protocol b)b)c) brings all trajectories to the globalminimum.[38] Eqs. (5-4) were Euler integrated with (cid:15) = 0 . M , γ l ( t =0) = − . M , ρ th = 0 . M, with D = 100 max(1 − ρ/ρ th , dt = 0 . M = max ≤ l ≤ N (cid:80) ¯Ω | A kl,i ,i , ··· ,i k | .The elements of dense (sparse) tensors were multipliedby 5 / N (15 / N ).[39] Eq. (6) was Euler integrated with dt = 0 .
001 with thesaem number of iterations as in other methods. The ini- tial conditions x l ( t = 0) were uniformly randomly dis-tributed in [ − . , . h j +1 = 1 . h j is updated eachtime (cid:80) l | dx l dt | < .
001 is satisfied.[40] G. Joya, M. Atencia, and F. Sandoval, “Hopfield neuralnetworks for optimization: study of the different dynam-ics,”
Neurocomputing , vol. 43, no. 1-4, pp. 219–237, 2002.[41] Eq. (7) was Euler integrated with dt = 0 .
02. iterationsequals 2000 timesteps. The initial conditions x l ( t = 0)are randomly distributed in [ − , β ( t = 0) = 10, while β j +1 = 0 . β j is updated each time (cid:80) l | dx l dt | < ..