Mean-field Coherent Ising Machines with artificial Zeeman terms

Coherent Ising Machine (CIM) is a network of optical parametric oscillators that solves combinatorial optimization problems by finding the ground state of an Ising Hamiltonian. In CIMs, a problem arises when attempting to realize the Zeeman term because of the mismatch in size between interaction and Zeeman terms due to the variable amplitude of the optical parametric oscillator pulses corresponding to spins. There have been three approaches proposed so far to address this problem for CIM, including the absolute mean amplitude method, the auxiliary spin method, and the chaotic amplitude control (CAC) method. This paper focuses on the efficient implementation of Zeeman terms within the mean-field CIM model, which is a physics-inspired heuristic solver without quantum noise. With the mean-field model, computation is easier than with more physically accurate models, which makes it suitable for implementation in field programmable gate arrays (FPGAs) and large-scale simulations. Firstly, we examined the performance of the mean-field CIM model for realizing the Zeeman term with the CAC method, as well as their performance when compared to a more physically accurate model. Next, we compared the CAC method to other Zeeman term realization techniques on the mean-field model and a more physically accurate model. In both models, the CAC method outperformed the other methods while retaining similar performance.


I. INTRODUCTION
The Ising model is useful in a wide variety of fields.Ising models, for instance, can not only demonstrate phase transitions of magnetic materials: ferro, ferri, anti-ferromagnets, and spin glasses [1][2][3][4][5][6] , but also capture other complex phenomena such as supercooled phases of water 7,8 , elections 9 , and disease spread 10,11 .Furthermore, Ising models can work as solvers for combinatorial optimization problems (COPs), such as quadratic unconstrained binary optimization, which are commonly encountered in real-life problems like scheduling problems 12,13 , portfolio optimization 14 , traffic volume optimization 15,16 , drug discovery 17 , and machine learning 18,19 .Therefore, it is possible to map many COPs to the problem of searching for the ground state of an Ising Hamiltonian.However, finding the ground state of an Ising Hamiltonian is a nondeterministic polynomial-time (NP-hard) problem, requiring computation time on an exponential scale relative to the size of the problem.There is a strong demand for Ising solvers for COPs due to the importance of these realworld problems.The development of dedicated hardware for searching the Ising Hamiltonian ground state has been active in recent years [20][21][22][23][24] .Particularly, Ising machines inspired by quantum physics have attracted wide attention due to their potential to overcome the above-mentioned difficulties for large COPs [25][26][27] .As of yet, however, there is still a problem with realizing Zeeman terms in quantum-inspired Ising machines.Ising Hamiltonian, including Zeeman term, is expressed as follows.
a) E. L. Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA b) School of Computing, Tokyo Institute of Technology, Tokyo, Japan where σ r represents the Ising spin variable taking either +1 or −1, and J rr ′ implies the coupling weight between spins r and r ′ .Here the diagonal entries are set to 0. The second term in eq.(1) indicates the external field present for each spin which is generally called the Zeeman term.Using the Ising Hamiltonian definition (eq.( 1)), the local field of each Ising spin can be described as follows.
Eq. ( 2) contains two terms, the first of which indicates the interaction term and the second of which indicates the Zeeman term.
Almost any COP can be characterized as an Ising problem with a Zeeman term 28 .Several such problems have been simulated with Ising solvers currently available, including Simulated Bifurcation Machine (SBM) in the traveling salesman problem 29 , Coherent Ising Machines (CIM) in the l 0 -regularization-based compressed sensing problem 27,30 , and Quantum Annealing in the Nursing Scheduling Problem 31 etc.Hence a Zeeman term must be implemented on quantuminspired Ising machines in order to apply them to solving realworld problems.However, due to the size mismatch between the interaction term and the Zeeman term, implementing a Zeeman term is difficult for machines with variable amplitude spins, such as CIMs and SBMs [32][33][34] .
In the conventional case the Ising Hamiltonian, which includes the Zeeman term (eq.( 1)), as well as the local field (eq.( 2)) for each Ising spin, is defined based on the assumption that each spin variable takes either +1 or −1.In CIMs, as well as SBMs, the spin amplitudes are continuous values that are different from ±1, so the interaction term and Zeeman term's contribution to the local field (e.g. for CIMs, eq. ( 6)) depend on the amplitude of the spin variable.As a consequence of the size mismatch between the interaction term and the Zeeman term, Ising machines may produce biased solutions which are inconsistent with the Hamiltonian defined by eq.(1).
For the search for ground states, CIM uses the minimumgain principle rather than thermal fluctuation as in classical annealing 35 and quantum fluctuation as in quantum annealing on D-wave and so on 36 .Based on a comparison between D-Wave (Chimera graph) and CIM (all-to-all coupling), it has been shown that CIMs demonstrate a performance advantage due to the differently implemented coupling 37 .On the other hand, there has been a claim that CIM has a better scaling capability when solving large-scale problems whereas SBM's performance is highly dependent on the hardware rather than the algorithm itself 38 .
The primary objective of this paper is the efficient implementation of Zeeman terms within mean-field CIM models that do not incorporate quantum noise terms and measurements (henceforth referred to as MFZ (Mean-Field-Zeeman)-CIM).The mean-field CIM model is a physics-inspired heuristic solver that does not accurately represent the CIM's behavior.However, due to their low computational costs, mean-field models are suitable for implementation with field programmable gate arrays (FPGAs) and for simulations on a large scale.So far, three approaches have been proposed to address the realization problem for CIM, namely the absolute mean amplitude method 33 , the auxiliary spin method 39 , and the chaotic amplitude control (CAC) method 40 .In this paper, we examine the applicability of CAC to realizing the Zeeman term in MFZ-CIM.Our results for the same optimization problem are compared to those of Ref. 40 's truncated-Wigner stochastic differential equations (SDEs) in order to demonstrate that MFZ-CIM performs almost as well as the truncated-Wigner SDEs.Using truncated-Wigner SDEs, we describe the approximate behavior of the experimental CIM device.Since truncated-Wigner SDEs exhibit the same performance in low quantum noise conditions as Positive-P SDEs, we focus only on this model 40 .Furthermore, this paper discusses the performance difference between these two models and the other Zeeman term realization methods as well.

A. Zeeman term implementations on CIM
Absolute Mean Amplitude: In the early stages of CIM research, it was suggested that the Zeeman term could be efficiently incorporated into the injection field by scaling it with the absolute mean of the amplitudes of the OPO pulses 33,41 .
Here, I in j,r is the injection field for the OPO pulse r. x r states the normalized in-phase amplitude of the OPO pulse r while J rr ′ and h r are the coupling weight and the Zeeman term described above.ζ is the adjustment parameter for the strength of the Zeeman term.Here j represents the feedback strength.
Auxiliary Spin: Recent efforts have been made to implement Zeeman terms in CIM and determine the ground state in an efficient manner 39,40,42 .It has been demonstrated by Singh et al., that the Zeeman terms in CIMs can be realized by using auxiliary spins 39 .In this case, the Zeeman term is incorporated into the two-body interaction term through the introduction of auxiliary spins to be included in a product with the Zeeman term as follows.
where x N+1 is an auxiliary amplitude to match the size of the Zeeman term to the interaction term and to transform the Zeeman term to an additional interaction term.As indicated in eq. ( 3), the injection field is reformulated only by the interaction term given in J rr ′ ∈ R (N+1)×(N+1) and x ∈ R (N+1) .The extended coupling matrix can be constructed by giving additional column and row vectors as J rN+1 = ζ h r and J N+1r ′ = ζ h r ′ , and taking J N+1,N+1 = 0. Currently, CIMs only support two-body interactions, which makes this method effective.Refer to Ref. 39 for a detailed explanation.
Chaotic Amplitude Control: In Ref. 40 , the Zeeman term has been implemented into the injection field through a technique known as Chaotic Amplitude Control (CAC).CAC is a technique that was proposed by Leleu et al., to overcome the problem of amplitude inhomogeneity in CIMs 43 .With CAC, the amplitudes of OPO pulses are forced to equalize to a set target value while forcefully correcting inhomogeneities resulting in a chaotic behavior which may result in escaping from local minima in the energy landscape 43 .By scaling the Zeeman terms with target amplitude to match the interaction term, Inui et al., in Ref. 40 proposed an efficient approach for implementing Zeeman terms in CIM as follows.
Here the target amplitude is indicated as τ.x is the in-phase amplitude of the OPO pulse, and e r is the auxiliary variable for the error feedback in the CAC feedback loop.Using two CIM models expressed as the Wigner stochastic differential equation (W-SDE) and the Positive-P stochastic differential equation (P-P-SDE), a high success probability of finding the ground state of a Sherrington-Kirkpatrick (SK) Hamiltonian with a random external field has been found as a result of implementing the Zeeman term with CAC 40 .

B. Mean-Field CIM model
There is one bottleneck in Ref. 40 regarding large-scale simulations, in that the more physically accurate SDEs used are quite difficult for a digital device such as an FPGA to implement.It is possible, however, to consider an alternative and a simpler differential equation (DE) model by disregarding the quantum noise present and any measurement effects.It is commonly referred to as the mean-field equations in literature 38,44,45 .The amplitude dynamics of mean-field DEs are governed by the equation below.
On the right-hand side of the equation, the first, second, and third terms represent linear loss, pump gain, and nonlinear saturation, respectively.As for the fourth term, it corresponds to the mutual coupling term.In comparison with eq. ( 7) and the amplitude governing SDE in Ref. 40 (Truncated Wigner and Positive-P) (see Appendices A and B), these relaxations result in a simplified deterministic CIM model which is more appropriate for use with digital hardware.Based on previous work from Ng et al., 45 , we introduced Gaussian white noise into the mean-field model at its initial state with a variance of 10 −4 while maintaining the advantages of the deterministic model.Utilizing the Mean-Field CIM (MF-CIM), this heuristic approach was used to construct a simple quantum-inspired Ising algorithm.In this case, the Wigner-type stochastic differential equation's (eq.(B4) -(B7) in Appendix B) amplitude variable is normalized with the saturation parameter and its quantum noise is ignored 45 .

III. NUMERICAL EXPERIMENTS
Numerical simulations were conducted in order to assess the effectiveness of CAC on MFZ-CIM.In this case, we are considering random SK Hamiltonians with a randomly generated Zeeman term.J rr ′ is constructed as a symmetric matrix where the elements are constructed using a random normal distribution with a mean of 0 and a variance of 1.The Zeeman term was also constructed by using a random normal distribution with a mean of 0 and a variance of 1.The diagonal axis of J rr ′ was set to 0. In order to compare with the solution energy produced by the CIM, the exact solution for each randomly generated J rr ′ and h r was calculated brute-force.Here we also consider Inui et al.,'s CIM model that is expressed as W-SDE and implemented Zeeman term with CAC (hereafter referred to as GATW (Gaussian-Approximation-Truncated-Wigner)-CIM) (given in eq. ( B4)-(B7) in Appendix B) and MFZ-CIM.Throughout all GATW-CIM simulations, the pump rate was scheduled according to the following schedule.
For each time-step t, p(t) indicates the calculated pump rate.As suggested in Ref. 40 , the target τ scheduling for GATW-CIM was as follows.
In this case, the calculated p(t) is used in order to calculate the τ(t) for each time-step.The saturation parameter g 2 corresponds to the quantum noise present in the GATW-CIM (see Appendices A and B).Throughout all MFZ-CIM simulations, the pump rate was constant as p = 0.57.Since MFZ-CIM does not take into account quantum noise present in the CIM, we use a simple linear τ(t) scheduling as follows.
In this case, τ(t) is linearly increased by each time-step t. τ 0 implies the starting target values while τ n represents the last target value after n time-steps.The final target value is given by (τ 0 + τ n ).We set τ 0 = 1 and τ n = 2. Ξ indicates the maximum time-steps assigned to the simulation.Throughout all simulations, feedback value e r is initially set to 1 in both GATW-CIM and MFZ-CIM.Consequently, the CACfeedback does not operate when β = 0. P sc indicates the average success probability in the figures.In a simulation, success is determined by the difference between the final energy calculated from the estimated final spin configuration by CIM and the brute-force calculation being less than 10 −4 .In order to calculate the Psc, we divided the number of successful runs by the number of simulations conducted overall.

A. Typical behavior of MFZ-CIM
An illustration of the typical amplitude x r evolution for an N=16 (here N refers to the system size) mean-field CIM with no Zeeman term (ζ = 0) and no CAC feedback (β = 0) can be found in Fig. 1a.In this case, the amplitudes are clearly in-homogeneous.In Fig. 1b, a random Zeeman term is introduced with ζ = 1 in an MFZ-CIM.There is, however, no CAC-feedback (β = 0) -hence, the amplitudes are still inhomogeneous.Fig. 1c depicts MFZ-CIM amplitudes with a random Zeeman term (ζ = 1) and a CAC feedback signal (β = 10).With CAC, evolving amplitudes have a chaotic nature until 4 photon life-times and after that amplitudes become homogeneous.An illustration of the corresponding error variables e r is provided in Fig. 1d.It is the abrupt jumps in error variables that cause chaotic fluctuations in amplitudes.

B. Relationship between performance and quantum noise
When the MFZ-CIM neglects quantum noise, we must ask whether this has an impact on performance for the better or for the worse as compared to the GATW-CIM.To understand this, we calculated the average success probability for an N=16 CIM-produced final spin states using 500 random SK Hamiltonians where the exact solution energy is calculated brute-force to compare.The red solid and dashed lines in Fig. 2 represent the results obtained by both CIMs without CAC-feedback (β = 0).The blue solid and dashed lines, on the other hand, represent the results with the CAC feedback (β = 10).In Fig. 2 the upper graph depicts the GATW-CIM results, while the lower graph illustrates the MFZ-CIM results.Zeeman term strength is increased on the horizontal axis, while success probability (P sc ) is indicated on the vertical axis.
It is evident that the performance of GATW-CIM increases when g 2 is small (g 2 = 10 −7 ) for β = 10 (blue solid line).In contrast, performance decreases when g 2 is increased to 10 −3 (blue dashed line).Comparing these two scenarios without CAC-feedback (β = 0), performance falls dramatically in both situations (red solid and dashed lines respectively).Meanwhile, in the MFZ-CIM case, the success probability is almost the same as that of the GATW-CIM for g 2 = 10 −7 (blue solid lines of Fig. 2a and Fig. 2b).Furthermore, MFZ-CIM tends to be slightly more successful in areas with a higher ζ value when β = 0 (red solid lines of Fig. 2a and Fig. 2b).
A further comparison is shown in Fig. 3 between the performance of GATW-CIM and MFZ-CIM when β is increased gradually with ζ .In this experiment, β is increased in intervals of 0.5 from 0.5 to 10. Figures 3a and 3b represent the GATW-CIM results and the MFZ-CIM results, respectively.The color plot of Fig. 3 indicates whether the success probability is higher or lower depending on β and ζ .Blueish hues indicate a lower success rate, while yellowish hues indicate a high success rate.Both MFZ-CIM and GATW-CIM exhibit similar performance, as can be seen from their respective figures.Even so, it is clear that MFZ-CIM tends to perform better than GATW-CIM in higher-ζ regions.For both models, j = 1 was given, while g 2 = 10 −7 used in GATW-CIM.On the horizontal axis, we increased ζ , and on the vertical axis, we indicated the probability of success.Using the solid blue, red and yellow lines, we can identify the CAC-feedback method in Fig 4 (eq.( 5) and eq.( 6)), absolute mean (eq.( 3)) and auxiliary spin (eq.( 4)) techniques, respectively.In both CIM models, the CAC-feedback method has a higher success probability compared to the other two methods.Note that both Absolute Mean Method and Auxiliary Spin are open-loop implementations.This means there is no dynam-ical modulation to the injection field I in j,r .CAC, however, is a closed-loop method in which error variable e r modulates I in j,r dynamically and individually.There is a significant decrease in the success probability of the auxiliary spin method in the higher ζ region, while both the CAC method and absolute mean method maintain relatively better performance.As a whole, both CIM models perform quite similarly for the respective Zeeman term realization methods.(a) Average P sc for GATW-CIM.(b) Average P sc for MFZ-CIM.The blue, red and yellow solid lines represent the CAC feedback, Absolute mean and Auxiliary spin method respectively.CAC strength was set to β = 10 with j = 1.Average success probabilities were calculated using 500 random SK Hamiltonians for N = 16 CIM models.

A. Effect of Zeeman term on performance
The acquired results suggest that both GATW-CIM and MFZ-CIM perform similarly with CAC-feedback present.
Without CAC, MFZ-CIM tends to perform better compared to when g 2 = 10 −3 for GATW-CIM.A possible explanation for this may be the quantum noise present in GATW-CIM.
Despite the mentioned similarities, there are also some slight differences.The performance of GATW-CIM is slightly superior in lower ζ situations, as shown in Fig. 3.While considering the overall performance in higher ζ cases, MFZ-CIM has a better performance.It is possible that these slight differences can be attributed to the effects of quantum noise in GATW-CIM.
Here because the interaction term becomes matched with the Zeeman term around ζ = 1, we expected MFZ-CIM and GATW-CIM with CAC to have the highest success probability around ζ = 1.However, as Fig. 3 shows, the success probability increases monotonically as ζ increases.However, without CAC, both models have their success probability maximized around ζ = 1, as shown in Fig. 3.There is a possibility that CAC may cause the destabilization of local minima and the development of dynamical complexity to differ depending on the value of ζ .Detailed analysis of these factors is needed in the future.

B. Evaluation with larger system sizes
The model we used in the manuscript is referred to as the MF-CIM in the field of nonlinear quantum optics.This heuristic approach has been used to construct a simple quantuminspired Ising algorithm using the MF-CIM.It can be derived by normalizing the Wigner-type stochastic differential equation's (eq.(B4) -(B7) in Appendix B) amplitude variable with the saturation parameter and ignoring its quantum noise 45 .It is called a mean-field model because ignoring quantum noise is justified by the mean-field approximation 45 .Generally, mean-field approximations require larger system sizes.However, our numerical experiments are intended to demonstrate the effectiveness of the quantum-inspired Ising algorithm using the MF-CIM, not its validity with the mean-field approximation for neglecting quantum noise.
It has been the practice in many previous studies to evaluate the performance of proposed algorithms for the SK Hamiltonians under small system sizes (around N = 2 − 16) since brute-force searches were needed for finding the SK Hamiltonians' ground state 40,45,46 .We therefore only ran numerical simulations for N = 16 in order to brute-force search the ground state, similar to other studies.Several more studies may be necessary on problems such as Wishart Planted Instances 47 and Compressed Sensing 27,30 , which allow larger systems (such as N = 200 and N = 4000) to be used because ground states can be determined using statistical mechanics and algorithmic hardness can be adjusted with planted ground states, allowing for the expansion of system sizes.

C. Effects of Quantum noise on CIM
In our view, the major difference between GATW-CIM and MFZ-CIM is attributable to the absence of quantum noise in MFZ-CIM.By normalizing the SDEs of GATW-CIM and considering that g 2 is 0, the differential equations (DEs) of MFZ-CIM can be derived (See Appendix C).In the CIM system, 1/g 2 represents the photon number in the system, where g 2 reflects the amount of quantum noise present in the system.In our simulations, g 2 is 10 −7 and 10 −3 for the GATW-CIM.Therefore, it is safe to say that GATW-CIM operates with a modicum of quantum noise.On the other hand, in MFZ-CIM, g 2 is zero.Thus, the performance discrepancies must be mainly due to the small quantum noise present in the GATW-CIM.

D. Performance difference of Ising machines on the SK model
In this paper, we evaluate the performance of the MF-CIM for fully connected SK models.Metastable states in the SK model exhibit ultrametric organization in phase space.Because of this complicated phase structure, searching for the ground state with CIM or other Ising machines can be difficult.D-Wave uses quantum annealing (QA) as its operating principle, and reaching the ground state is accomplished by reducing the transverse field sufficiently slowly over time 48 .
In contrast, CIMs search for a ground state in a quantum parallel manner prior to reaching a threshold pump power, and once the threshold is reached, the ground state's amplitude is amplified according to the minimum gain principle 46 .In the context of fully connected SK models, CIM outperformed D-Wave experimentally 37 .Unlike D-Wave, which has an exponential computation time proportional to exp(O(N)), CIM has an exponential computation time proportional to exp(O( √ N)), where N is the size of the problem 38 .Furthermore, Kako et al (2020) demonstrated efficient sampling of ground states using truncated Wigner approximated SDEs with CAC on a fully connected Max-Cut model with N = 16 and 8, where the Max-Cut model can be defined by the same Hamiltonian as the SK model 46 .This paper shows that under implementing the Zeeman term, MF-CIM with CAC achieves almost the same performance as truncated Wigner SDEs with CAC on success probability.

E. Future work on MFZ-CIM
It has been demonstrated by Gunathilaka et al., that NP-Hard problems such as l 0 -regularized compressed sensing (L0RCS) can be solved with CIMs that employ target amplitude to implement Zeeman terms.However, the SDEs employed in Ref. 27 are difficult to implement in digital hardware, such as FPGAs.With the simplicity of MFZ-CIM, we plan to explore numerically the possibility of implementing it directly on digital hardware and applying it to L0RCS.

VI. CONCLUSION
In this paper, we have evaluated the mean-field model of CIM with a Zeeman term present.It is apparent from our acquired results that the performance of the MFZ-CIM with CAC-feedback is roughly similar to the more accurate SDEs to the physical CIM called GATW-CIM with CAC-feedback.Even though it was relatively hard to introduce the Zeeman term to the CIM because of the mismatch in size between interaction and Zeeman terms, the introduction of CAC has enabled a way to realize the Zeeman terms while keeping relatively good performance.According to our results, the CAC technique is more effective than the previous Zeeman term realization technique, for the mean-field CIM model.Furthermore, the use of simplified SDEs with a Zeeman term will be beneficial when implementing CIM on digital hardware for solving real-world problems.
A CIM consists of a set of optical parametric oscillators (OPOs), where oscillations above the threshold limit constitute an optimum solution to a Hamiltonian 35 .Using an optical cavity in conjunction with a phase-sensitive amplifier (PSA), a coherent state can be achieved in the CIM.This allows the spin-up state to be defined as 0-phase and the spin-down state as π-phase, and the overall state to be defined as an Ising spin model.The J matrix or the mutual coupling matrix is created using a mutual injection field.The measurement feedback CIM (MFB-CIM) master equation can be expressed as follows 40 .

FIG. 4 .
FIG. 4. Performance difference with other Zeeman term realization techniques.(a)Average P sc for GATW-CIM.(b) Average P sc for MFZ-CIM.The blue, red and yellow solid lines represent the CAC feedback, Absolute mean and Auxiliary spin method respectively.CAC strength was set to β = 10 with j = 1.Average success probabilities were calculated using 500 random SK Hamiltonians for N = 16 CIM models.