## Abstract

As a solver for non-deterministic polynomial time (NP)-hard combinatorial optimization problems, the coherent Ising machine (CIM) is in the early stages of research, and the potential of this innovative physical system will be developed. Here, we propose a speed-up coherent Ising machine with a squeezed feedback system, which we call S-CIM. We couple squeezed feedback pulses generated by the squeezed feedback system into the degenerate optical parametric oscillator (DOPO) network. Simulations indicate that quantum inseparability of the coupled DOPO network is further enhanced during the whole optimization process, and quantum fluctuations are significantly smaller around the oscillation threshold. Computation experiments are performed on MAX-CUT problems of order between 4 and 20000. Numerical results demonstrate that S-CIM increases the optimal normalized output by 2.27% and significantly reduces the optimal computation time by 75.12%.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Many combinatorial optimization problems in computer science [1], drug research [2], protein folding [3], and big data processing [4] belong to the non-deterministic polynomial time (NP)-hard class [5]. The traditional von-Neumann computer has no enough capability to solve these problems. It can't find the exact solution with satisfying efﬁciency. Recently, some interesting attempts [6–9], such as simulated annealing generated by thermal annealing program in simulated metallurgical process and quantum annealing based on quantum tunneling theory, are proposed to tackle the NP-hard problems approximately. However, how to control the tight connection between quantum bits more efficiently needs to be explored [10]. These methods have inspired us to find alternatives from other reliable physical schemes.

Recently, Stanford University proposed and implemented a coherent Ising machine (CIM) based on degenerate optical parametric oscillators (DOPOs) [11,12] for seeking the optimal solutions of the NP-hard problems. The basic operating mechanism of CIM is derived from mapping the global photon decay rate of the DOPO network to the Ising Hamiltonian [5]

where*J*indicates the coupling coeﬃcient between the

_{jl}*j*and the

_{th}*l*spins, and

_{th}*σ*= ±1 denotes two spin states, |↑> and |↓>. Due to the inherent characteristics of CIM, this experimental setup can implement an artiﬁcial Ising spin system, where the Ising couplings are realized by injecting feedback pulses into the network of DOPO. Furthermore, Ising spin states are suitably represented with binary phase states, |0 > and |π>, of the above-threshold DOPO pulses. When the gain of CIM gradually approaches the global minimum loss by increasing the external pump, the DOPO network will begin to oscillate and the phase state (spin state) configuration will be in the exact or approximate ground state of the Ising Hamiltonian. Finding the ground state of the Ising Hamiltonian is classified as the Ising problem, which also belongs to the NP-hard class [13]. CIM is used as an efficient optimization solver because of the reducibility of the NP problems [14]. Recently, Stanford University and NTT improved the scheme of the CIM by adding quantum measurement feedback control system, and successfully reduced the size of the machine [15]. Shortly afterward, NTT used dual-pump four-wave mixing in a highly nonlinear fiber [16] to generate more than 50,000 time-multiplexed DOPOs [17], which indicates that the DOPO network can be used as a stable large-scale artificial spin system. However, there are still some factors that reduce the optimization efficiency of the CIM. The signal-to-noise ratio (SNR) of the measurement by balanced homodyne detector was disturbed due to the antisqueezing effect generated in the fiber ring cavity [18] and the vacuum fluctuations introduced by the beam splitter [15]. In addition, the turn-on-delay oscillation effect postpones the end of the optimization process [19].

_{j }In this paper, a speed-up coherent Ising machine with a squeezed feedback system (S-CIM) is proposed. We formulate a calculation model of the coupled DOPO network based on squeezed feedback process by using truncated Wigner representation. Then we obtain the c-number stochastic differential equations (CSDEs) for the network consisting of *N* DOPOs. To make the results clear, a simple model composed of two coupled DOPOs is taken to study the quantum properties of the S-CIM. When the squeezed feedback pulses are injected into the 2-DOPO network, the variances of the quadrature amplitude component for DOPOs decreases, and quantum inseparability is enhanced throughout the optimization process. The simulations about the evolutions of the spin state (phase state) configurations and the average photon number in the fiber ring cavity have consistent results, suggesting that S-CIM oscillates earlier than CIM to get optimization results. Finally, the computational performance of the S-CIM is numerically evaluated by solving large-scale MAX-CUT problems of order up to 20000.

## 2. Calculation model

The schematic of the S-CIM shown in Fig. 1(a) consists of three parts: a *N*-DOPO network as an artiﬁcial Ising spin system; a squeezed feedback system consisting of a subthreshold optical parametric oscillator, a phase sensitive amplify (PSA), a balanced homodyne detector, an analog-to-digital converter (ADC), a ﬁeld-programmable gate array (FPGA), a digital-to-analog converter (DAC), an intensity modulator (IM) and a phase modulator (PM); and an injection coupler for realizing Ising couplings.

A portion of the S-CIM external pump is injected in the fiber ring cavity to form *N* signal DOPO pulses, and another portion is pumped into the subthreshold optical parametric oscillator to generate squeezed vacuum states. In the squeezed feedback process, the *j _{th}* signal fields are extracted from the beam splitter 1 (BS

_{1}) to measure its quadrature amplitude component ${\tilde{x}_j}$. The output of the balanced homodyne detector is converted into a digital signal and input into the FPGA. Then we use an electronic digital circuit to compute the

*l*

_{th}feedback pulse corresponding to the

*j*

_{th}DOPO pulse. The output signal drives the IM and PM to accurately achieve the

*l*

_{th}feedback pulse $\mathop \sum \nolimits_l {\xi _{jl}}{\tilde{x}_l}$, where ${\xi _{jl}}$ is proportional to the coupling coefficient

*J*

_{jl}, and

*J*is a $N \times N$ symmetric matrix, with

_{jl}*J*= 0. Afterward the squeezed vacuum state is injected into the feedback pulse, the squeezed feedback pulse with a given Ising coupling

_{jj}*J*emerge. Finally, the squeezed feedback pulses are coupled into the ring cavity through the beam splitter 2 (BS

_{jl}_{2}), which produces a coupled DOPO network. Each part of the calculation model will be analyzed in detail, and the truncated Wigner representation will be used to eventually obtain the c-number stochastic differential equations (CSDEs) of the S-CIM containing

*N*-coupled DOPOs.

#### 2.1 Degenerate optical parametric oscillators

The S-CIM is driven by a master laser at frequency *ω _{s}*. The external pump

*F*with a frequency of

_{p}*ω*is excited by the second harmonic generation (SHG). The external pump

_{p }= 2ω_{s}*F*is divided into two optical paths via a beam splitter. One enters the fiber ring cavity and the other injects the subthreshold optical parametric oscillator. Due to the second-order nonlinear effect of the crystal, the pump entering the ring cavity generates signal pulses sequence at frequency

_{p}*ω*.

_{s}*N*DOPOs run along the fiber ring cavity at regular time intervals

*t*′ =

*T*′/

*N*, where

*T*′ is the roundtrip time when a DOPO pulse propagates through the cavity. Therefore, the number of independent DOPO can be added by increasing the repetition frequency of the pump pulse or by extending the length of the ring cavity. Assuming that only two DOPOs are coupled in the cavity, the total Hamiltonian [11,20] of the network is obtained:

*H*

_{free}describes the free field Hamiltonian of signal fields and pump fields in the 2-DOPO network inside the cavity, where $\hat{a}_{sj}^{\dagger} $, ${\hat{a}_{sj}}$ are the creation operator and annihilation operator of the

*j*signal field and $\hat{a}_{pj}^{\dagger} $, ${\hat{a}_{pj}}$ are the correspondences for the pump field.

_{th}*H*

_{int}represents interaction Hamiltonian, where

*к*is the parametric gain derived from a second-order nonlinear crystal.

*H*

_{pump}denotes pump field Hamiltonian excited by the external pump fields, where $\epsilon $ is pump rate. At last,

*H*

_{sr}shows the interaction between reservoir and the cavity ﬁelds. It consists of two signal and two pump fields, where ${\gamma _p}$, ${\gamma _s}$ are the pump and signal optical decay rates and ${\hat{\Gamma }}_{Rsj}^{\dagger} $, ${\hat{\Gamma }}_{Rpj}^{\dagger} $ are reservoir operators.

#### 2.2 Squeezed feedback pulse

A small portion of the *j _{th}* signal DOPO pulse is extracted from the fiber ring cavity by the BS

_{1}. Its quadrature amplitude component is measured by the balanced homodyne detector. The transmittance rate of BS

_{1}is defined as

*T*

_{1}= sin

^{2}

*θ*

_{1}. Since the measured intensity (1-

*T*

_{1}) is relatively small, the quantum coherence between different amplitudes can be improved [18]. When the pulse propagates through the BS

_{1}, it combines with the vacuum fluctuations of the external interference, resulting in measurement errors of the detector. The annihilation operator of the vacuum fluctuation is defined as a

_{vac}, and the output coupled field operator is

*ξ*[11] is achieved.

Then the measurement result of homodyne detector is transmitted to FPGA, we can calculate the feedback pulse for ${\hat{a}_{sj}}$. After the calculation, the output electrical signal drives the intensity and phase modulator to modulate the extracted pulse ${\hat{a}_{sj,out}}$ to feedback pulse ${{\xi }_{jl}}{\hat{a}_{sl}}$.

At the same time, the external pump generated by the master laser propagates into the subthreshold optical parametric oscillator to produce squeezed vacuum states [23,24]. After the squeezed vacuum state is injected [25,26] into the measurement feedback circuit, the obtained annihilation operator of the squeezed feedback pulse satisfies

*r*is squeeze parameter. Since the evolution of the network only depends on the quadrature amplitude components of the signal DOPO pulses, we suppress the quantum noise of the quadrature amplitude component at the cost of the increased the quantum noise of the quadrature phase component.

#### 2.3 Squeezed feedback process and the coupled DOPO network

In the squeezed feedback process, the signal DOPO pulse and squeezed feedback pulse, the former of which is prepared in a coherent state, are combined by the BS_{2.} The transmittance rate of the BS_{2}, defined as *T*_{2} = sin^{2}*θ*_{2}, is set very high (*T*_{2} ≈ 1). In this parameter region, the quantum noise is nearly completely suppressed because of the injected squeezed vacuum state. Carrying out this coupling procedure requires three valid modes [27] of the pulses as shown in Fig. 1(b): ${\hat{a}_{sj}}$ (Black arrow) is Mode 1^{+}, ${\hat{b}_{sl}}$ (yellow arrow) is Mode 2^{+}, and ${\hat{c}_{sjl}}$ (red arrow) is Mode 1^{−}. Mode 1^{+} describes the DOPOs in the main cavity which are not extracted by the BS_{1}. Mode 2^{+} represents squeezed feedback pulses which are modulated by the squeezed feedback system. We assume that the state of the Mode 1^{−} corresponding to the coupled DOPO pulse is |*ψ*> = *S*_{2}(*ξ2*)*D*_{1}(*α _{sj}*)

*S*

_{1}(

*r*

_{1})|0>, where

*α*is the complex amplitude of the jth signal DOPO pulse. ${D_1}$ is the displacement operator of Mode 1

_{sj}^{+},

*S*

_{1}is the squeezed operator of Mode 1

^{+}, and

*S*

_{2}is the squeezed operator of Mode 2

^{+}. Here, let

*r*

_{1}= 0, so Mode1

^{+}and Mode 2

^{+}are coherent state and squeezed state, respectively. Then, the creation and annihilation operator of Mode1

^{−}are obtained as follows:

*W*(

*α*) [29] for the four modes to expand the field density operator

*ρ*:

*A*

_{s}_{12}=

*gα*

_{s}_{12}is the normalized complex amplitude, $g = \frac{\kappa }{{\sqrt {2{\gamma _s}{\gamma _p}} }}$ is the saturation parameter, $E = \frac{\kappa }{{{\gamma _s}{\gamma _p}}}\epsilon $ is the normalized pump rate,

*ξ*is the effective coupling coefficient from

_{jl}*l*squeezed feedback pulse to

_{th}*j*DOPO pulse,

_{th}*τ*=

*γ*is the normalized time, and the noise term is

_{s}t*N*-coupled DOPO network:

*m*=

*j*and

*n*=

*l*cannot be established at the same time to avoid double calculation of the coupled pulse.

## 3. Quantum properties

The quantum inseparability and fluctuations of 2-coupled DOPO network are discussed based on the truncated Wigner function [29] because of its convenience of evaluating the expectation value of a symmetrically ordered operator:

*α*} = (

*α*

_{s}_{12},

*α*

_{s}_{21})

*for the 2-coupled DOPO network. In the present case, the quadrature amplitude component and the quadrature phase component of the coupled DOPO pulse are ${\hat{x}_{sjl}} = \frac{{({\hat{c}_{jl}^{\dagger} + {{\hat{c}}_{jl}}} )}}{2}$ and ${\hat{p}_{sjl}} = \frac{{({{{\hat{c}}_{jl}} - \hat{c}_{jl}^{\dagger} } )}}{{2i}}$, respectively. The EPR-type operators ${\hat{u}_ + } = {\hat{x}_{s12}} + {\hat{x}_{s21}}$ and ${\hat{v}_ - } = {\hat{p}_{s12}} + {\hat{p}_{s21}}$ are used to evaluate quantum inseparability. The evaluation criteria of quantum inseparability is ${\Delta }\hat{u}_ + ^2 + {\Delta }\hat{v}_ - ^2 < 1$ [30]. As shown in Fig. 2(a), the total variances of the EPR-type operators for the CIM with rapidly increasing pump rate, the CIM with gradually increasing pump rate and the S-CIM with gradually increasing pump rate are compared which are calculated by the truncated Wigner representation. The pumping schedule of the rapidly increasing pump rate is that the pump rate linearly increases from zero to 1.5 times DOPO oscillation threshold within the normalized time*

^{T}*τ*= 200, i.e.

*E*= 1.5(

*τ*/200), and the pumping schedule of the gradually increasing pump rate is

*E*= 1.5(

*τ*/800). Here the oscillation threshold refers to the threshold pump rate $E_{th}^{(0 )} = 1$ of a solitary DOPO. As shown in Fig. 2(a), the CIM with rapidly increasing pump rate has a weak quantum inseparability below the coupled oscillation threshold

*E*

_{th}= 1-ξ [11]. Here, we set the coupling constant ξ = 0.6 and the saturation parameter is

*g*= 0.01 for better results. When the coupled DOPO network evolves above the coupled oscillation threshold, the quantum inseparability gradually disappears, i.e. ${\Delta }\hat{u}_ + ^2 + {\Delta }\hat{v}_ - ^2 \approx 1$. There is a small peak from

*E*= 0.54 to

*E*= 0.72 due to the turn-on-delay oscillation effect, which leads to the reduction of the optimization efficiency. Note that the small peak disappears in the CIM with gradually increasing pump rate, which indicates that the slower pumping schedule can weaken the turn-on-delay oscillation effect. As for the S- CIM with gradually increasing pump rate, it can be seen that the turn-on-delay oscillation effect is weakened, and the quantum inseparability is further deepened over the entire pumping schedule. Moreover, we compare the quantum inseparability of coupled DOPO network based on two different feedback process models: the valid coupling modes in this paper (Fig. 2(a)) and the unitary displacement operator in the Heisenberg picture $D(\beta {\theta _2}) = exp(\beta {\theta _2}{a^ + } - {\beta^\ast }{\theta _2}a)$ [18] (Fig. 2(b)), where $|\beta > $ is the feedback pulse prepared in a squeezed state. The numerical values confirm that the difference in the total variances $\left\langle {\Delta \widehat u_ +^2} \right\rangle + \left\langle {\Delta \widehat v_ -^2} \right\rangle $ computed by two methods is within the statistical error.

The variances of the quadrature amplitude component represent quantum fluctuations which include the anti-squeezed effect caused by PSA inside the cavity and the vacuum fluctuations introduced from the BS_{1}. Large fluctuations will reduce the SNR of detection, cause errors in the squeezed feedback system and affect the accuracy of the final optimization result. As shown in Fig. 3(a), the variances of the quadrature amplitude component for CIM and S-CIM are compared under the slowly pumping schedule *E* = 1.5(*τ*/800). Below the coupled oscillation threshold, the variances begin to increase from the vacuum state (0.5) with the rise of the pump rate. Once the coupled network oscillates, the variance gradually decreases to the original value. Comparing with CIM, S-CIM further reduces the quantum fluctuations in the coupled DOPO network. Similarly, we compared the quantum fluctuations of the coupled DOPO network obtained by the two feedback process models. As shown in Figs. 3(a) and 3(b), the difference in the variances of the quadrature amplitude component ${\widehat x_{s12}}$ computed by two methods is still within the statistical error. Therefore, the simulation results obtained by the valid coupling modes also have acceptable accuracy.

## 4. Optimization process

Through the analysis of the optimization output statistics, the optimization process of the CIM is divided into four stages [19]: quantum parallel search, quantum ﬁltering, spontaneous symmetry breaking, and quantum-to-classical crossover. Then in the 2-coupled DOPO network, the four specific optimization stages of the S-CIM are compared with the CIM.

The spin states (phase states) of two coupled DOPOs are mapped to four possible spin configurations |↑↓>, |↑↑>, |↓↓> and |↓↑>. During the optimization process of S-CIM, when the coupled DOPO network oscillates, the ground state of the Ising Hamiltonian is found, and the spin configuration is the solution of the corresponding Ising problem. Figures 4(a)–4(d) present the post-selected probabilities of four spin states during the S-CIM evolution as the final selected configuration is |↑↓>. As Fig. 4(a) shows, two DOPOs are independent of each other, and four spin configurations have the same probability of 0.25 which implies that the S-CIM system is at an early stage of ‘quantum parallel search’. Next, as shown in Fig. 4(b), the probability amplitudes of the spin configurations corresponding to the ground state, |↑↓> and |↓↑>, are amplified, and the probability amplitudes of the spin configurations corresponding to the excited state, |↑↑> and |↓↓>, are deamplified. The network evolves from linear superposition to quantum inseparability, which is called ‘quantum filtering’. Then, as shown in Fig. 4(c), the S-CIM system has a 50% probability of selecting one of the two ground states to amplify its probability amplitude. Here, the selected state is |↑↓>. Meanwhile, system deamplifies the probability amplitude of another unselected state |↓↑>. This process is called ‘spontaneous symmetry breaking’. Figure 4(d) shows the last step ‘quantum-classical crossover’. The probability amplitude of the selected state |↑↓> is increased to the maximum, while the probabilities of other spin states are reduced to the minimum. The optimization process ends, and the final optimization outputs are obtained by using S-CIM.

The probability evolutions of the two ground-state |↑↓> and |↓↑> are used to characterize the specific optimization stage, where the final selected state is |↑↓> and the unselected state is |↓↑>. Here, the CIM and the S-CIM have the same pumping schedule *E* = 1.5(*τ*/200). Firstly, let us focus on the CIM. When *τ* = 0, the two ground states |↑↓> and |↓↑> have the same probability of 0.25, corresponding to ‘quantum parallel search’ (green dashed line in Fig. 5). When 3 ≤ *τ* < 28, the probabilities of the two ground states gradually increases because of ‘quantum filtering’ (red dashed line in Fig. 5). When 28 ≤ *τ* < 84, the probability of the selected state |↑↓> continues to increase, but the probability of the unselected state |↓↑> starts to decrease, which are caused by ‘spontaneous symmetry breaking’ (blue dashed line in Fig. 5). Finally, the vivid result emerges at *τ* = 84, which indicates the stage of ‘quantum-classical crossover’. After calculation, the coupled oscillation threshold of the CIM is *E*_{th}* _{ }*= 0.4, and the corresponding normalization time is

*τ*

_{s}= 53. However, the CIM actually oscillates at

*τ*

_{d}= 84, which is due to the turn-on-delay oscillation effect. It can be clearly seen that in the three stages of ‘quantum filtering’, ‘spontaneous symmetry breaking’ and ‘quantum-classical crossover’, the time consumption of the S-CIM is significantly shorter than that of the CIM. The coupled DOPO network starts to oscillate at

*τ*= 71.

The dynamic oscillation thresholds *E*_{d} is accurately defined as the value of *E* that maximize $dlog(n )/dlog(E )$, where *n* is the average photon number in the fiber ring cavity. Figure 6 illustrates that S-CIM oscillates earlier than CIM according to the evolution of the average photon number $\langle{n}\rangle$, which is consistent with the analysis of the ground-state spin configurations (Fig. 5).

## 5. Computational experiments

The success probability and the computation time of CIM and S-CIM in solving the NP-hard MAX-CUT problems are used to evaluate their performance. Firstly, the simplest MAX-CUT-3 problem is calculated, which consists of four vertices and six edges with identical weight ${J_{ij}} ={-} 1$. The reason for using a small instance is that the accuracy of the device is verified by exhausting all possible outcomes. The optimized outputs of CIM and S-CIM are shown in Fig. 7. When the network has four coupled DOPOs, there will be sixteen possible outputs. Considering the degeneracy of the phase states (spin states), three of the eight spin configurations correspond to the ground state of the Ising Hamiltonian, {|↑↓↓↑>,|↑↑↓↓>,|↑↓↑↓>}, i.e. the solution of the simplest MAX-CUT-3 problem. In order to make the comparison between CIM and S-CIM more obvious, the feedback rate (The rate at which squeezed feedback pulse and signal DOPO pulse are coupled) is reduced, and the maximum success rate does not reach 100%. Figure 7 shows that S-CIM can increase the success rate of computation by 5.4% in the 4-DOPO network.

Next, large-scale instances with 800-20000 vertices and 0.02%-6% edge density are explored. The high order MAX-CUT problems were obtained from G-set graphs, which were randomly constructed by using a machine-independent graph generator created by G. Rinaldi [31]. The normalized optimization results and the average computation time are indicators for evaluating the performance of S-CIM. The computation time is defined as the evolution time of the quadrature amplitude component for DOPO from the initial state to the steady state.

When solving the same sample G1, Fig. 8 shows the evolution time of the quadrature amplitude component under two different schemes in one trial of numerical simulations. The bifurcation of quadrature amplitude components for CIM at 45 ≤ *τ* < 396 is caused by the spontaneous symmetry breaking. In the S-CIM, the stage of the spontaneous symmetry breaking occupies a shorter time range at 38 ≤ *τ* < 92. Finally, CIM evolves to steady state at *τ* = 396 [in Fig. 8(a)], while S-CIM reaches steady state at *τ* = 92 [in Fig. 8(b)]. This means that S-CIM completes the optimization earlier than CIM, which is in good agreement with the theoretical analysis.

Finally, the performance evaluation indicators of S-CIM are recorded in Table 1. *V* denotes the number of vertices in sample G-set graph, *E* denotes the number of edges, *C*_{GW} denotes the optimal output of the Goemans-Williamson (GW) algorithm based on semi-definite programming (SDP) [32], *U*_{SDP} denotes the solution to the SDP upper bound of the MAX-CUT problem, *C*_{CIM} and < *C*_{CIM}> denote the best and average normalized outputs of running the CIM 500 times, respectively. Similarly, *C*_{S-CIM} and < *C*_{S-CIM}> represent the best and average normalized outputs of running the S-CIM 500 times, respectively. In order to compare these schemes, every output is normalized by (*C *+ *E*_{neg})/(*U*_{SDP}+*E*_{neg}), where *E*_{neg} denotes the number of negative edges in sample G-set graph, *C* denotes the cut value of three schemes. <*τ*_{CIM}> and <*τ _{S-}*

_{CIM}> are the average computation time of running CIM and S-CIM 500 times, respectively. The average normalized solutions of S-CIM are about 0.18%-2.27% higher than that of CIM. Because the vacuum fluctuations in the ring cavity of S-CIM are smaller than that of CIM, the detection efficiency is improved, and the optimization outputs are more satisfactory. The reason why the improvement is not significant is that, as the number of vertices increases, the number of suboptimal solutions increases at a faster rate, greatly reducing the accuracy of the final results. In addition, the average computation time of CIM and S-CIM is compared in Table 1. For all large-scale instances, the computation time of S-CIM is 17.93%-75.12% less than that of CIM. Quantum inseparability accelerates the coupled DOPO network oscillating, which improves the optimization efficiency of the S-CIM system.

## 6. Conclusion

Squeezed feedback pulses are combined with signal DOPO pulses to generate an *N*-coupled DOPO network, which is used as a speed-up coherent Ising machine to solve the NP-hard MAX-CUT problem with shorter computation time and better optimization results. The calculation model of the S-CIM is formulated by using truncated Wigner representation, and its quantum properties are simulated. There is a turn-on-delay oscillation effect in the CIM. It leads to a difference between the dynamic threshold and the static threshold of the system. The coupled DOPO network actually takes more time to start oscillating. However, in the S-CIM with the schedule of the gradually increasing pump rate, the squeezing parameter of squeezed vacuum injection increases with increasing pump rate. Then the small peak indicating the turn-on-delay oscillation effect disappears, and the quantum inseparability of the system is enhanced. As a result, the time required for the network to search for the ground-state spin configuration is shortened, and the evolution of the average photon number in the fiber ring cavity is accelerated. All these suggest that S-CIM can get the optimization outputs earlier. In addition, the squeezed vacuum injection deamplifies the noise of the quadrature amplitude component at the cost of amplifying the noise of the quadrature phase component, which weakens the quantum fluctuations of the coupled DOPO network. Computation experiments of the large-scale MAX-CUT problems (*N *= 800-20000) demonstrate that the normalized optimization result for S-CIM is 2.27% higher than that of CIM in the best case, and the computation time is 75.12% less than that of CIM. These numerical results are consistent with the findings of the quantum properties, indicating that the introduction of the squeezed feedback system improves the optimization efficiency of CIM and gives a great push to the application of the CIM.

## Funding

National Natural Science Foundation of China (11604377, 61775234); Qingdao National Laboratory for Marine Science and Technology (QNLM2016ORP0111); CAS “Light of West China” Program (XAB2017B18); Natural Science Basic Research plan in Shaanxi Province of China (2018JQ6067); The Independent Project of State Key Laboratory of Transient Optics and Photonics.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **F. Barahona, “On the computational complexity of Ising spin glass models,” J. Phys. A: Math. Gen. **15**(10), 3241–3253 (1982). [CrossRef]

**2. **M. Marc, G. Parisi, and M. Virasoro, * Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications* (World Scientific Publishing Company, 1987).

**3. **D. B. Kitchen, H. Decornez, J. R. Furr, and J. Bajorath, “Docking and scoring in virtual screening for drug discovery: methods and applications,” Nat. Rev. Drug Discovery **3**(11), 935–949 (2004). [CrossRef]

**4. **J. D. Bryngelson and P. G. Wolynes, “Spin glasses and the statistical mechanics of protein folding,” Proc. Natl. Acad. Sci. **84**(21), 7524–7528 (1987). [CrossRef]

**5. **I. H. Witten, E. Frank, M. A. Hall, and C. J. Pal, * Data Mining: Practical Machine Learning Tools and Techniques* (Morgan Kaufmann, 2016).

**6. **S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science **220**(4598), 671–680 (1983). [CrossRef]

**7. **Z. Su, D. Xue, and Z. Ji, “Designing LED array for uniform illumination distribution by simulated annealing algorithm,” Opt. Express **20**(S6), A843–855 (2012). [CrossRef]

**8. **P. Honzatko, J. Kanka, and B. Vrany, “Retrieval of the pulse amplitude and phase from cross-phase modulation spectrograms using the simulated annealing method,” Opt. Express **12**(24), 6046–6052 (2004). [CrossRef]

**9. **T. Kadowaki and H. Nishimori, “Quantum annealing in the transverse Ising model,” Phys. Rev. E **58**(5), 5355–5363 (1998). [CrossRef]

**10. **E. G. Rieffel, D. Venturelli, B. O’Gorman, M. B. Do, E. M. Prystay, and V. N. Smelyanskiy, “A case study in programming a quantum annealer for hard operational planning problems,” Quantum Inf. Process. **14**(1), 1–36 (2015). [CrossRef]

**11. **Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, “Coherent Ising machine based on degenerate optical parametric oscillators,” Phys. Rev. A **88**(6), 063853 (2013). [CrossRef]

**12. **A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, “Network of time-multiplexed optical parametric oscillators as a coherent Ising machine,” Nat. Photonics **8**(12), 937–942 (2014). [CrossRef]

**13. **M. Gu and Á Perales, “Encoding universal computation in the ground states of Ising lattices,” Phys. Rev. E **86**(1), 011116 (2012). [CrossRef]

**14. **M. R. Garey and D. S. Johnson, * Computers and Intractability: a Guide to the Theory of NP-Completeness* (Freeman, 2009).

**15. **Y. Haribara, S. Utsunomiya, and Y. Yamamoto, “Computational principle and performance evaluation of coherent Ising machine based on degenerate optical parametric oscillator network,” Entropy **18**(4), 151–166 (2016). [CrossRef]

**16. **T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, “Large-scale Ising spin network based on degenerate optical parametric oscillators,” Nat. Photonics **10**(6), 415–419 (2016). [CrossRef]

**17. **H. Takesue and T. Inagaki, “10 GHz clock time-multiplexed degenerate optical parametric oscillators for a photonic Ising spin network,” Opt. Lett. **41**(18), 4273–4276 (2016). [CrossRef]

**18. **A. Yamamura, K. Aihara, and Y. Yamamoto, “Quantum model for coherent Ising machines: Discrete-time measurement feedback formulation,” Phys. Rev. A **96**(5), 053834 (2017). [CrossRef]

**19. **D. Maruo, S. Utsunomiya, and Y. Yamamoto, “Truncated Wigner theory of coherent Ising machines based on degenerate optical parametric oscillator network,” Phys. Scr. **91**(8), 083010 (2016). [CrossRef]

**20. **K. Takata, A. Marandi, and Y. Yamamoto, “Quantum correlation in degenerate optical parametric oscillators with mutual injections,” Phys. Rev. A **92**(4), 043821 (2015). [CrossRef]

**21. **D. F. Walls, “Squeezed states of light,” Nature **306**(5939), 141–146 (1983). [CrossRef]

**22. **C. M. Caves, “Quantum limits on noise in linear amplifiers,” Phys. Rev. D **26**(8), 1817–1839 (1982). [CrossRef]

**23. **H. P. Yuen, “Two-photon coherent states of the radiation field,” Phys. Rev. A **13**(6), 2226–2243 (1976). [CrossRef]

**24. **G. J. Milburn and D. F. Walls, “Squeezed states and intensity fluctuations in degenerate parametric oscillation,” Phys. Rev. A **27**(1), 392–394 (1983). [CrossRef]

**25. **L. A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, “Generation of squeezed states by parametric down conversion,” Phys. Rev. Lett. **57**(20), 2520–2523 (1986). [CrossRef]

**26. **Z. Dutton, J. H. Shapiro, and S. Guha, “LADAR resolution improvement using receivers enhanced with squeezed-vacuum injection and phase-sensitive amplification,” J. Opt. Soc. Am. B **27**(6), A63–A72 (2010). [CrossRef]

**27. **C. M. Caves, “Quantum-mechanical noise in an interferometer,” Phys. Rev. D **23**(8), 1693–1708 (1981). [CrossRef]

**28. **H. J. Carmichael, * Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations* (Springer-Verlag, 2002).

**29. **D. F. Walls and G. J. Milburn, * Quantum Optics* (Springer, 2007).

**30. **L. M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, “Inseparability criterion for continuous variable systems,” Phys. Rev. Lett. **84**(12), 2722–2725 (2000). [CrossRef]

**31. **C. Helmberg and F. A. Rendl, “spectral bundle method for semidefinite programming,” SIAM J. Control **10**(3), 673–696 (2000). [CrossRef]

**32. **M. X. Goemans and D. P. Williamson, “Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming,” J. Assoc. Comput. Mach. **42**(6), 1115–1145 (1995). [CrossRef]