## Taking the pulse of optimization

Finding the optimum solution of multiparameter or multifunctional problems is important across many disciplines, but it can be computationally intensive. Many such problems defined as computationally difficult can be mathematically mapped onto the so-called Ising problem, which looks at finding the minimum energy configuration for an array of coupled spins. Inagaki *et al.* and McMahon *et al.* show that an optical processing approach based on a network of coupled optical pulses in a ring fiber can be used to model and optimize large-scale Ising systems. Such a scalable architecture could help to optimize solutions to a wide range of complex problems.

## Abstract

The analysis and optimization of complex systems can be reduced to mathematical problems collectively known as combinatorial optimization. Many such problems can be mapped onto ground-state search problems of the Ising model, and various artificial spin systems are now emerging as promising approaches. However, physical Ising machines have suffered from limited numbers of spin-spin couplings because of implementations based on localized spins, resulting in severe scalability problems. We report a 2000-spin network with all-to-all spin-spin couplings. Using a measurement and feedback scheme, we coupled time-multiplexed degenerate optical parametric oscillators to implement maximum cut problems on arbitrary graph topologies with up to 2000 nodes. Our coherent Ising machine outperformed simulated annealing in terms of accuracy and computation time for a 2000-node complete graph.

Combinatorial optimization is an important task for various applications such as drug discovery, resource optimization in wireless networks, and machine learning (*1*). Many such problems are known to be nondeterministic polynomial time (NP)–hard or NP-complete problems, which are considered difficult to solve efficiently with modern digital computers. Several heuristic or “greedy” algorithms have been proposed in the field of computer science to find the optimal or near-optimal solutions to those problems. Among such algorithms, simulated annealing (SA) has performed well in solving various types of combinatorial optimization problems (*2*). It is known that combinatorial optimization problems can be mapped onto ground-state search problems of the Ising model with polynomial resources (*3*). Several approaches have been proposed and demonstrated to find solutions to Ising problems using networks of artificial spins based on various physical systems (*4*–*7*).

A coherent Ising machine (CIM) is one such artificial Ising machine, in which a laser (*8*) or a degenerate optical parametric oscillator (DOPO) (*9*) represents an artificial spin. This proposal was motivated by early studies on the quantum behaviors of a DOPO (*10*–*14*), and proof-of-principle CIM experiments have already been reported using time-multiplexed DOPOs (*15*, *16*). To date, physical Ising machines (*4*–*7*) have suffered from a limited number of spin-spin couplings due to localized spin implementations. For example, in a recent quantum annealing system based on 1152 superconducting qubits, the number of spin-spin couplings was 3360, constrained by the chimera graph structure (*17*). This limited number of spin-spin couplings has made it difficult to map various real-world optimization problems to the hardware. Here, we report a CIM that consists of 2048 DOPOs with full spin-spin couplings. A measurement and feedback (MFB) scheme enabled us to implement all possible connections among 2048 spins, which amounts to 2,096,128 spin-spin couplings, making it possible to solve arbitrary graph problems.

An Ising Hamiltonian without an external magnetic field is given by (1)(*18*), where σ* _{i}* ∈ {–1,+1} denotes the

*i*th Ising spin,

*J*is the interaction strength between the

_{ij}*i*th and

*j*th spins, and

*N*is the total number of spins. The Ising machines described above are designed to find the ground state or low-energy states of the Ising Hamiltonian. A DOPO takes only a 0 or π phase relative to the pump phase above the threshold, so it can be used as a stable artificial spin realized with photonics technologies (

*19*–

*22*). A CIM solves the Ising problem according to the minimum-gain principle (

*8*). According to this principle, if the optical couplings between the DOPOs are adjusted to implement {

*J*}, the DOPO network oscillates in the phase configuration with the lowest loss, which corresponds to the ground state of the Ising Hamiltonian in Eq. 1.

_{ij}The size of the problem that a CIM can handle is determined by the number of DOPOs and their couplings. Delayed interferometers were used to couple the DOPOs in the previous CIM experiments (*15*, *16*); however, this scheme does not scale as the problems become more complex. To avoid this issue, Haribara *et al*. proposed an MFB scheme that makes it possible to couple any combination of *N* DOPOs (*23*).

Figure 1 shows a schematic of a CIM with MFB (*24*). The CIM uses a periodically poled lithium niobate (PPLN) waveguide, which is placed in a 1-km fiber cavity, as a phase-sensitive amplifier (PSA) that amplifies only the 0 or π phase components relative to the pump phase (*25*, *26*). Other components in the 1-km fiber cavity include a 9:1 coupler (coupler 1) for extracting the DOPO pulses for balanced homodyne detection (BHD), a 1-km dispersion-shifted fiber, and another 9:1 coupler (coupler 2) for injecting feedback signal pulses. When we start pumping the PSA, quadrature-squeezed noise pulses are generated as a result of spontaneous parametric downconversion. The noise pulses undergo repeated phase-sensitive amplification, which leads to the formation of DOPO pulses. Because the cavity round-trip time is 5 μs and the pump pulse interval is 1 ns, the cavity can accommodate a maximum of 5082 DOPOs, of which 2048 DOPOs are used as Ising spins (henceforth “signal DOPOs”). Another 2834 DOPOs, which we call “training DOPOs,” are always being oscillated and are used to obtain error signals for phase-locking the 1-km fiber cavity, the local oscillator for the BHD and the coupling pulses (*24*). As a result, we can perform sequential computations stably for a few seconds. The remaining 200 DOPOs are turned off (i.e., the pump pulses are turned off at the corresponding time slots) to separate the signal and training DOPO pulses.

To perform a computation, we turn the signal DOPOs on and off by controlling the power of the pump pulses for the signal DOPOs in a 5-ms period (fig. S2A). During the computation, the in-phase components of the signal DOPOs from coupler 1 are measured using the BHD for every circulation of the DOPOs in the cavity. The power of the local oscillator for the BHD is 2.5 mW, with which the shot noise power is ~8 dB larger than the thermal noise power, and thus a shot noise–dominant coherent measurement is achieved. The measured is supplied to a field-programmable gate array (FPGA) module. The given problem, specified by the spin-spin interaction terms {*J _{ij}*} (a 2048 × 2048 adjacency matrix), is input in advance into the FPGA module. The current FPGA module only accepts

*J*values of {–1, 0, 1}. With and {

_{ij}*J*}, the FPGA module calculates the feedback signal to the

_{ij}*i*th pulse

*f*= for every DOPO circulation, where

_{i}*r*denotes a constant that determines the coupling strength and was set at ~0.01 in the following experiments. A push-pull modulator imposes the feedback signal

*f*on the coupling pulses, which are extracted from the pump pulses. Then, the modulated coupling pulses are launched into the cavity through coupler 2 so that a coupling pulse that conveys

_{i}*f*is injected into the

_{i}*i*th DOPO. Note that this MFB sequence is repeated throughout the entire computational sequence. With the current FPGA module, it takes about ~60 s to transfer the problem matrix and the computation result data between the FPGA module and the computer, because of the use of a slow serial communication interface (RS-232C) that is the only available port in the module. The data transfer time can be much shorter than the computation trial time if we make use of high-speed interfaces (

*24*).

We studied the computational accuracy of the CIM and compared it with SA run on a single core of an Intel Xeon X5650 (2.67 GHz) processor for 2000-node graphs. With the CIM and SA, we solved maximum-cut (MAX-CUT) problems (*24*) for three types of 2000-node undirected graphs: G22, G39, and K_{2000}. G22 and G39 are sparsely connected graphs taken from the G-set (*27*). G22 is a random graph with a 1% edge density, and G39 is a scale-free graph in which the node degree distribution follows power-law scaling. K_{2000} is a fully connected complete graph with 1,999,000 undirected edges, which is randomly weighted by {–1, +1}. These three instances were randomly constructed with a machine-independent graph generator (*28*). We also solved the same graphs using the semidefinite programming relaxation algorithm proposed by Goemans and Williamson (GW-SDP), because this algorithm guarantees to deliver solutions to MAX-CUT problems with an expected value of at least 87.8% of the optimal value (*29*). We used the cut values obtained by GW-SDP as an index for the computational accuracy.

An example of the temporal evolution of the signal DOPO amplitudes observed while solving the K_{2000} graph is shown in Fig. 2A, which portrays the cooperative behavior of the DOPOs leading to the final stable states after frequent changes of their signs at the beginning of the computation. To show the detailed dynamics of each signal DOPO, we plot the measured in-phase components of 50 of the 2000 signal DOPOs in Fig. 2B. As shown in the inset of Fig. 2B, we observed the evolution of the DOPO amplitudes searching for a lower-energy state almost immediately after the pumping started; the dynamics was complex, especially in the first 200 circulations. This result confirmed a reduction in the Ising energy (fig. S4).

Using the CIM, we performed 100 sequential computations for MAX-CUT on the three graphs. The cut values obtained with the CIM and SA for each graph are summarized in Fig. 3. Here, the final DOPO states for each computation trial were stably recorded sequentially 100 times, thanks to the phase stabilization of the optical system. The SA result was obtained in 50 ms for each trial without the data read and write time. With GW-SDP, a computation trial took about 1 s. For all three graphs, the CIM achieved larger cut values than those obtained by GW-SDP, which confirmed that a physical system based on the CIM concept successfully worked to find solutions to the large-scale Ising problems with as many as 2000 spins. SA produced better values than the CIM when solving sparse graphs G22 and G39, as confirmed by the cut histograms in Fig. 3. In contrast, the cut obtained by the CIM for the K_{2000} graph exceeded that obtained by SA with one-tenth of the computation time needed for SA.

We then compared the computation time needed to find an approximate solution for the MAX-CUT of the K_{2000} graph. Because the exact solution of this problem is not known, we cannot evaluate the computation time to find the optimal solution. Hence, we measured the computation times needed with the CIM and SA to reach the target Ising energy obtained with GW-SDP, namely –60,278 (or 29,619 cuts). The temporal evolution curves for the CIM and SA are shown in Fig. 4 by blue and black curves, respectively. Here, we performed the computation 100 times; the curves indicating the shortest times needed to reach the target energy are shown for both the CIM and SA. The times needed to reach the target Ising energy were about 70 μs for the CIM and 3.2 ms for SA, implying that the CIM could obtain the benchmark Ising energy more quickly than with SA by a factor of nearly 50. This result suggests that a CIM can outperform SA on a CPU when solving dense graph problems.

According to the numerical study of Maruo *et al*. (*14*), the formation of a superposed state is observed when an optically coupled CIM is operated in the shot noise–limited regime. This suggests that the CIM may exhibit quantum behavior, although experimental confirmation is yet to be achieved. A more detailed discussion on the quantumness of the CIM is provided in (*24*).

## Supplementary Materials

www.sciencemag.org/content/354/6312/603/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S5

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵See supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**Supported by the Impulsing Paradigm Change Through Disruptive Technologies (ImPACT) Program of the Council of Science, Technology and Innovation (Cabinet Office, Government of Japan). We thank H. Nishimori for fruitful discussions, K. Inaba for fruitful discussions, and H. Tamura for various types of support during this research. S.U. and H. Takesue are inventors on patent application PCT/JP2015/059057 submitted by the National Institute of Informatics (NII) and Nippon Telegraph and Telephone (NTT) Corporation that covers the coherent Ising machine based on the measurement and feedback scheme. A.M. and S.U. are inventors on patent application PCT/US2014/046025 submitted by Stanford University and NII that covers the implementation of a coherent Ising machine using degenerate optical parametric oscillators. T.U., K.E., and H. Takenouchi are inventors on patent application PCT/JP1012/000360 submitted by NTT that covers the configurations of phase-sensitive amplifiers based on periodically poled lithium niobate waveguides.