Speed-up coherent Ising machine with a squeezed feedback system

Lihaonan Luo; Lihaonan Luo; Hongjun Liu; Hongjun Liu; Nan Huang; Zhaolu Wang

doi:10.1364/OE.381850

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]

(1)$$H ={-} \mathop \sum \limits_{1 \le j < l \le N} {J_{jl}}{\sigma _j}{\sigma _l},$$

where J_jl indicates the coupling coeﬃcient between the j_th and the l_th spins, and σ_j= ±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].

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.

Fig. 1. (a). The schematic of the speed-up coherent Ising machine with a squeezed feedback system (S-CIM). The PSA connected to BS₁ amplifies the quadrature amplitude component of each extracted DOPO pulse, while the BS₂ couples the modulated squeezed feedback pulse with the target DOPO pulse to implement the given Ising Hamiltonian. A: the network of N degenerate optical parametric oscillators (DOPOs). B: the squeezed feedback system. C: the injection coupler. (b) The corresponding calculation model of the coupled DOPO network based on squeezed feedback process.

Download Full Size | PDF

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₁) 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_jl is a $N \times N$ symmetric matrix, with J_jj = 0. Afterward the squeezed vacuum state is injected into the feedback pulse, the squeezed feedback pulse with a given Ising coupling J_jl emerge. Finally, the squeezed feedback pulses are coupled into the ring cavity through the beam splitter 2 (BS₂), 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_p with a frequency of ω_p= 2ω_s is excited by the second harmonic generation (SHG). The external pump F_p 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 ω_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:

(2)$$H = {H_{\textrm{free}}} + {H_{\textrm{int}}} + {H_{\textrm{pump}}} + {H_{\textrm{sr}}},$$

(3)$${H_{\textrm{free}}} = \hbar {\omega _s}\mathop \sum \limits_{j = 1}^2 \hat{a}_{sj}^{\dagger} {\hat{a}_{sj}} + \hbar {\omega _s}\mathop \sum \limits_{j = 1}^2 \hat{a}_{pj}^{\dagger} {\hat{a}_{pj}},$$

(4)$${H_{\textrm{int}}} = \frac{{i\hbar \kappa }}{2}\mathop \sum \limits_{j = 1}^2 ({\hat{a}{{_{sj}^{\dagger} }^2}{{\hat{a}}_{pj}} - \hat{a}{{_{pj}^{\dagger} }^2}{{\hat{a}}_{sj}}} ),$$

(5)$${H_{\textrm{pump}}} = i\hbar \sqrt {{\gamma _p}} \mathop \sum \limits_{j = 1}^2 ({\epsilon \hat{a}_{pj}^{\dagger} {e^{ - i{\omega_p}t}} - \epsilon {{\hat{a}}_{pj}}{e^{i{\omega_p}t}}} ),$$

(6)$${H_{\textrm{sr}}} = i\hbar \mathop \sum \limits_{j = 1}^2 {\hat{a}_{sj}}{\hat{\Gamma }}_{Rsj}^{\dagger} \sqrt {{\gamma _s}} + {{\hat{\Gamma }}_{Rsj}}\hat{a}_{sj}^{\dagger} \sqrt {{\gamma _s}} + {\hat{a}_{pj}}{\hat{\Gamma }}_{Rpj}^{\dagger} \sqrt {{\gamma _p}} + {{\hat{\Gamma }}_{Rpj}}\hat{a}_{pj}^{\dagger} \sqrt {{\gamma _p}} ,$$

Here, 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_th signal field and $\hat{a}_{pj}^{\dagger} $, ${\hat{a}_{pj}}$ are the correspondences for the pump field. 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₁. Its quadrature amplitude component is measured by the balanced homodyne detector. The transmittance rate of BS₁ is defined as T₁ = sin²θ₁. Since the measured intensity (1- T₁) is relatively small, the quantum coherence between different amplitudes can be improved [18]. When the pulse propagates through the BS₁, 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

(7)$$\hat{a}_{sj}^{\prime} = \sqrt {{T_1}} {\hat{a}_{sj}} - \sqrt {1 - {T_1}} {\hat{a}_{vac}},$$

The PSA [21,22] in front of the balanced homodyne detector is used to amplify the quadrature amplitude components of the out-coupled field:

(8)$${\hat{a}_{sj,out}} = G\left( {\sqrt {{T_1}} {{\hat{a}}_{sj}} - \sqrt {1 - {T_1}} {{\hat{a}}_{vac}}} \right) + \sqrt {G - 1} \left( {\sqrt {{T_1}} \hat{a}_{sj}^ +{-} \sqrt {1 - {T_1}} \hat{a}_{vac}^ + } \right),$$

where G is the gain of the PSA. the SNR of the measurement is improved and a large coupling constant ξ [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

(9)$${\xi _{jl}}{\hat{b}_{sl}} = \frac{{{\xi _{jl}}}}{{G\sqrt {{T_1}} }}({{{\hat{a}}_{sj,out}}\cosh (r )- \hat{a}_{sj,out}^{\dagger} {e^{i\theta }}\sinh (r )} ),$$

where 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₂, defined as T₂ = sin²θ₂, is set very high (T₂ ≈ 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₁. 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)D₁(α_sj)S₁(r₁)|0>, where α_sj is the complex amplitude of the jth signal DOPO pulse. ${D_1}$ is the displacement operator of Mode 1⁺, S₁ is the squeezed operator of Mode 1⁺, and S₂ is the squeezed operator of Mode 2⁺. Here, let r₁ = 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:

(10)$$\begin{array}{l} \hat{c}_{sjl}^{\dagger} = \sqrt {{T_2}} \hat{b}_{sl}^{\dagger} + \sqrt {1 - {T_2}} \hat{a}_{sj}^{\dagger} \\ {{\hat{c}}_{sjl}} = \sqrt {{T_2}} {{\hat{b}}_{sl}} + \sqrt {1 - {T_2}} {{\hat{a}}_{sj}}, \end{array}$$

In this 2-DOPO network, the master equation for the fields consisting of two signal modes and two pump modes can be obtained by the standard technique [28]. Then, we use the Wigner function W(α) [29] for the four modes to expand the field density operator ρ:

(11)$$\rho = \int {{e^{{\lambda ^{\ast }}\hat{c} - \lambda \hat{c}}}^{^{\dagger} }} \left\{ {\int {{e^{\lambda {\alpha^{\ast }} - {\lambda^{\ast }}\alpha }}} W(\alpha )d\alpha } \right\}d\lambda ,$$

where $\hat{c} = {({{{\hat{c}}_{s12}},{{\hat{c}}_{s21}},{{\hat{a}}_{p12}},{{\hat{a}}_{p21}}} )^T}$ and ${\lambda} = ({{\lambda_{s12}},{\lambda_{s21}},{\lambda_{p12}},{\lambda_{p21}}} )$. Substituting Eq. (11) into the master equation to obtain the Fokker–Planck equation, and the third and higher-order terms of the Fokker–Planck equations are truncated to obtain a set of CSDEs by the Ito rules [29]:

(12)$$\begin{array}{l} \textrm{d}{\alpha _{s12}} = ({ - {\gamma_s}{\alpha_{s12}} + \kappa {\alpha_{p12}}\alpha_{s12}^{\ast }} )dt + \sqrt {{\gamma _s}} d{W_{s12}}(t )\\ \textrm{d}{\alpha _{s21}} = ({ - {\gamma_s}{\alpha_{s21}} + \kappa {\alpha_{p21}}\alpha_{s21}^{\ast }} )dt + \sqrt {{\gamma _s}} d{W_{s21}}(t )\\ \textrm{d}{\alpha _{p12}} = \left( { - {\gamma_p}{\alpha_{p12}} - \frac{\kappa }{2}\alpha_{s12}^2 + \epsilon } \right)dt + \sqrt {{\gamma _p}} d{W_{p12}}(t )\\ \textrm{d}{\alpha _{p21}} = \left( { - {\gamma_p}{\alpha_{p21}} - \frac{\kappa }{2}\alpha_{s21}^2 + \epsilon } \right)dt + \sqrt {{\gamma _p}} d{W_{p21}}(t ), \end{array}$$

where ${W_X}(t )$ is the Wiener increment. Under the condition of ${\gamma _p} \gg {\gamma _s}$ [11], the pump modes are adiabatically eliminated. Then, the CSDEs are reduced to normalized signal stochastic differential equations:

(13)$$\begin{array}{l} \textrm{d}{A_{s12}} = \{{ - {A_{s12}} + ({E - A_{s12}^2} )A_{s12}^{\ast} - {\xi_{21}}{A_{s21}}} \}d\tau + gdW_{s12}^{\prime}(\tau )\\ \textrm{d}{A_{s21}} = \{{ - {A_{s21}} + ({E - A_{s21}^2} )A_{s21}^{\ast} - {\xi_{12}}{A_{s12}}} \}d\tau + gdW_{s21}^{\prime}(\tau ), \end{array}$$

where A_s₁₂ = gα_s₁₂ 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, ξ_jl is the effective coupling coefficient from l_th squeezed feedback pulse to j_th DOPO pulse, τ = γ_st is the normalized time, and the noise term is

(14)$$\begin{array}{l} dW_{s12}^{\prime} = \frac{1}{4}{e^{ - \textrm{E}{\gamma _S}}}\textrm{d}{W_{s12}}(\tau )+ {A_{s12}}\textrm{d}{W_{p1}}(\tau )\\ dW_{s21}^{\prime} = \frac{1}{4}{e^{ - \textrm{E}{\gamma _S}}}\textrm{d}{W_{s21}}(\tau )+ {A_{s21}}\textrm{d}{W_{p2}}(\tau ), \end{array}$$

Finally, Eqs. (13) and (14) can be extended to an N-coupled DOPO network:

(15)$$\begin{array}{l} d{A_{sjl}} = \left\{ { - {A_{sjl}} + ({E - A_{sjl}^2} )A_{sjl}^{\ast} - \mathop \sum \limits_{m = 1}^N \mathop \sum \limits_{n = 1}^N {\xi_{mn}}{A_{smn}}} \right\}d\tau + gdW_{sjl}^{\prime}(\tau )\\ dW_{sjl}^{\prime}(\tau )= \frac{1}{4}{e^{ - \textrm{E}{\gamma _S}}}\textrm{d}{W_{sjl}}(\tau )+ {A_{sjl}}\textrm{d}{W_{pj}}(\tau )({j \ne l,m \ne n} ). \end{array}$$

where 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:

(16)$${\left\langle {\hat{c}_{s12}^{{\dagger} j}\hat{c}_{s21}^{{\dagger} k}\hat{c}_{s12}^l\hat{c}_{s21}^m} \right\rangle _s} = \int {\alpha _{s12}^{{\ast}j}\alpha _{s21}^{{\ast }k}\alpha _{s12}^l\alpha _{s21}^m} W({\{\alpha \}} )d\alpha ,$$

here, {α} = (α_s₁₂, α_s₂₁)^T 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 τ = 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.

Fig. 2. Total variances of the EPR operators ${\Delta }\hat{u}_ + ^2 + {\Delta }\hat{v}_ - ^2$ versus normalized pump rate E for three different CIM schemes. The numerical values are computed from two feedback process models: the valid coupling modes approach (Fig. 2(a)) and the unitary displacement operator approach (Fig. 2(b)).

Download Full Size | PDF

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₁. 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.

Fig. 3. Variances of the quadrature amplitude component ${\hat{x}_{s12}}$ versus normalized pump rate E for CIM and S-CIM with gradually increasing pump rate, respectively. The numerical values are computed from two feedback process models: the valid coupling modes approach (Fig. 3(a)) and the unitary displacement operator approach (Fig. 3(b)). Since the curves of ${\hat{x}_{s12}}$ and ${\hat{x}_{s21}}$ are almost identical, only the former is given here.

Download Full Size | PDF

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.

Fig. 4. In the 2-coupled DOPO network, the probabilities of finding |↑↓>, |↑↑>, |↓↓> and |↓↑> states at four different stages of the optimization process when the final selected ground state is |↑↓>. The gain of the coupled DOPO network (yellow dotted line) gradually increases from below the coupled oscillation threshold to reach the minimum global loss (blue curve) corresponding to the ground state of the Ising Hamiltonian.

Download Full Size | PDF

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.

Fig. 5. In the 2-coupled DOPO network, the probability evolutions of two ground-state spin configurations versus normalized time τ for CIM and S-CIM. The green, red, blue and pink portions in each curve represent the four procedures of the optimization process of quantum parallel search, quantum filtering, spontaneous symmetry breaking, and quantum-to-classical crossover, respectively. The optimization process of CIM is divided by yellow dotted lines, and the optimization process of S-CIM is divided by yellow solid lines.

Download Full Size | PDF

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).

Fig. 6. In the 2-coupled DOPO network, the average photon number n versus normalized pump rate $p = \frac{E}{{{E_d}}}$ for CIM and S-CIM, respectively.

Download Full Size | PDF

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.

Fig. 7. In the 4-coupled DOPO network, the comparison of spin configuration distributions in 1000 trials of numerical simulations for CIM and S-CIM, respectively.

Download Full Size | PDF

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.

Fig. 8. Quadrature amplitude components of DOPOs versus normalized time τ for CIM and S-CIM against the MAX-CUT problem on the G1 graph.

Download Full Size | PDF

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.

Table 1. Performance comparison of CIM and S-CIM in computing large-scale MAX-CUT problems on sample G-set graphs.

View Table

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]

Graph	V	E	U_SDP	C_GW	<C_CIM>	C_CIM	<C_S-CIM>	C_S-CIM	<τ_CIM>	<τ_S-CIM>
G1	800	19176	12083	94.57%	95.23%	95.92%	95.67%	96.04%	406	101
G6	800	19176	2656	94.48%	94.91%	95.23%	95.22%	95.78%	382	114
G11	800	1600	629	93.27%	94.15%	94.46%	95.45%	95.67%	379	107
G18	800	4694	1166	92.82%	95.22%	95.40%	95.51%	95.72%	412	138
G22	2000	19990	14136	91.91%	94.88%	94.92%	95.50%	95.61%	544	201
G27	2000	19990	4141	91.74%	93.16%	95.21%	95.43%	95.74%	571	263
G35	2000	11778	8014	92.93%	95.20%	95.39%	95.38%	95.54%	526	557
G39	2000	11778	2877	92.27%	94.08%	94.64%	95.21%	95.39%	489	308
G43	1000	9990	7032	92.92%	94.12%	94.52%	95.15%	95.66%	476	285
G51	1000	5909	4006	93.34%	93.28%	94.06%	94.48%	94.82%	436	314
G55	5000	12948	11039	90.06%	93.69%	94.37%	94.73%	95.21%	762	399
G57	5000	1000	3885	92.37%	94.72%	94.93%	94.98%	95.02%	536	376
G60	7000	17148	15222	89.89%	94.16%	95.07%	94.81%	95.13%	783	404
G64	7000	41459	10466	91.43%	94.38%	94.78%	94.89%	95.25%	574	421
G67	10000	20000	7744	92.15%	93.54%	93.89%	94.62%	94.93%	697	450
G70	10000	9999	9683	96.33%	93.47%	94.11%	94.47%	94.76%	685	513
G81	20000	40000	15656	91.95%	93.21%	93.85%	94.36%	94.41%	803	659

Speed-up coherent Ising machine with a squeezed feedback system

Abstract

1. Introduction

2. Calculation model

2.1 Degenerate optical parametric oscillators

2.2 Squeezed feedback pulse

2.3 Squeezed feedback process and the coupled DOPO network

3. Quantum properties

4. Optimization process

5. Computational experiments

6. Conclusion

Funding

Disclosures

References

Cited By

Figures (8)

Tables (1)

Equations (16)

Optics Express