## 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

Unconventional, special-purpose machines may aid in accelerating the solution of some of the hardest problems in computing, such as large-scale combinatorial optimizations, by exploiting different operating mechanisms than those of standard digital computers. We present a scalable optical processor with electronic feedback that can be realized at large scale with room-temperature technology. Our prototype machine is able to find exact solutions of, or sample good approximate solutions to, a variety of hard instances of Ising problems with up to 100 spins and 10,000 spin-spin connections.

Combinatorial optimization problems, including many nondeterministic polynomial-time–hard (NP-hard) problems, are central in numerous important application areas, including operations and scheduling, drug discovery, finance, circuit design, sensing, and manufacturing. Despite large advances in both algorithms and digital computer technology, even typical instances of NP-hard problems that arise in practice may be very difficult to solve on conventional computers. There is a long history of attempts to find alternatives to current von Neumann–computer–based methods for solving such problems, including use of neural networks realized with analog electronic circuits (*1*, *2*) and by using molecular computing (*3*–*5*). Both lines of investigation continue to inspire related work (*6*, *7*). A major topic of contemporary interest is the study of adiabatic quantum computation (AQC) (*8*) and quantum annealing (QA) (*9*, *10*). Sophisticated AQC/QA devices are already under study (*11*–*14*), but providing dense connectivity between qubits remains a major challenge (*15*), with important implications for the efficiency of AQC/QA systems (*16*).

Networks of coupled optical parametric oscillators (OPOs) are an alternative physical system, with an unconventional operating mechanism (*17*–*20*), for solving the Ising problem (*21*, *22*) and by extension many other combinatorial optimization problems (*23*). Formally, the *N*-spin Ising problem is to find the configuration of spins (*i* = 1, ..., *N*) that minimizes the energy function , where the particular problem instance being solved is specified by the *N* × *N* matrix *J* (with elements *J _{ij}*) and the length-

*N*vector

*h*(with elements

*h*).

_{i}We have realized a system with a scalable architecture that uses measurement feedback in place of a network of optical delay lines [which were used in initial, low-connectivity, nonreprogrammable demonstrations of the concept (*18*, *24*, *25*)]. Our 100-spin Ising machine allows connections between any spin and any other spin and is fully programmable. We show that measurement-feedback–based OPO Ising machines can solve many different Ising problems, and in cases in which exact solutions are not easy to obtain, we can find good approximate solutions.

The schematic of our experimental setup (Fig. 1) shows that our Ising machine is formed by the combination of time-division–multiplexed OPOs (*18*) in a single fiber-ring cavity, with measurement and feedback (injection) stages that act to couple the pulses in the cavity such that the Ising Hamiltonian is realized. Details are provided in the supplementary materials (*26*).

We have focused our experiments on Ising problems on undirected, unweighted graphs (*V*, vertices; *E*, edges), where *J _{ij}* = –1 when spin

*i*and spin

*j*are connected, and

*J*= 0 otherwise, and for which the linear (Zeeman) terms

_{ij}*h*are zero. An Ising problem of this form is equivalent to the Max-Cut problem on the underlying graph (

_{i}*17*). Max-Cut is the problem of partitioning the vertices of a graph into two disjoint subsets such that the number of edges between the two subsets is maximized; the partition is called a cut. Max-Cut remains NP-hard even when the input is restricted to unweighted cubic graphs (

*27*). We refer interchangeably to the Max-Cut problem and to the Ising problem, and the solutions thereof. The energy of a particular spin configuration {σ

*}*

_{i}

_{i}_{=1,...,}

*for an Ising problem is given by . The same spin configuration {σ*

_{N}*} represents a cut of size ; there is a direct mapping between the energies in the Ising problem and the cut sizes in the Max-Cut problem, and minimizing the Ising energy maximizes the cut.*

_{i}The unweighted, undirected Möbius ladder (cubic) graph with *N* = 16 vertices is shown in Fig. 2A. We programmed the corresponding *J* matrix into our feedback electronics and ran the problem on our experimental apparatus. Figure 2B shows, from a single run of this single problem instance, the evolution of the in-phase component of each OPO pulse (*c _{i}*) as a function of the number of times each pulse circulated around the cavity (the number of round trips). The computation time is given by

*T*

_{comp}=

*T*

_{rt}

*N*

_{rt}, where

*T*

_{rt}= 1.6 μs is the single round trip time and

*N*

_{rt}is the number of round trips. Each spin σ

*in the Ising problem is represented by a single OPO pulse; if the in-phase component*

_{i}*c*of the

_{i}*i*th OPO pulse is less than zero (

*c*< 0), we make the spin assignment σ

_{i}*= –1, and if the in-phase component is greater than zero (*

_{i}*c*> 0), we make the spin assignment σ

_{i}*= +1. As the feedback signal level is gradually increased, the OPO amplitudes increase. Most of the OPOs obtain their ultimate signs by round trip 60. By round trip 100, all the OPOs have reached a steady state. Shown in Fig. 2C is the graph cut size, or equivalently the Ising energy, represented by the spin configuration after each roundtrip for the run shown in Fig. 2B. We see that the system evolves toward solutions representing larger graph cuts, or equivalently lower Ising energies, and the steady-state configuration is that of a ground state. The system finds the ground state within 120 μs. If we rerun the computation many times, we find that our system always returns the optimal solution for this particular problem instance (Fig. 2D). The distributions of obtained solutions for two other*

_{i}*N*= 16 cubic graphs are shown in Fig. 2, E and F. In both cases, the system finds the ground state in a large fraction of the runs, and when the system does not return a ground state, it returns a state that has energy close to that of the ground state, a good approximate solution. To illustrate that these instances are not special cases for which the Ising machine finds exact solutions (ground states), we attempted to solve all possible problem instances of

*N*= 16 cubic graphs, of which there are 4060. We were able to find ground states with probability greater than 20% for every single instance (Fig. 2G). In this experiment, and all that follow, every run was set to proceed for

*N*

_{rt}= 300 round trips, a computation time of

*T*

_{comp}= 480 μs.

We next explored how the Ising machine performs as a function of the size of the problem. In this first scaling study, we used the Möbius ladder graphs with varying size. This is a convenient choice because there is a closed-form solution for the maximum cut and Ising energy for every instance in this family of graphs. The Möbius ladder graphs with *N* = 8, 12, ..., 100 were solved (Fig. 3A). For each instance, we performed multiple sets of 100 runs. Each 100-run set resulted in a solution-energy histogram (an example for *N* = 8 is shown in Fig. 3B). For the runs in which the ground state is not found, the system again finds good approximate solutions (Fig. 3, C and D).

To understand how the performance of our system scales as a function of problem size for more general problems, 10 random instances of cubic graphs were generated for each graph size *N* = 10, 20, ..., 100, and each instance run on our apparatus. The success probabilities for each instance are shown in Fig. 4A, aggregated by graph size. In contrast to the previous results, not only the probability of obtaining a ground state but also the probability of obtaining a solution that is within *x*% of the optimal (maximum cut) solution is shown. The error bars indicate a spread in success probabilities primarily due to the fact that different random instances of each size may be more or less difficult to solve.

Figure 4B shows, for each graph size, the total computation time required to obtain a solution with the given accuracy, with probability 99% within that time, and is derived directly from Fig. 4A. The total time is given by *T*_{comp} × ⌈log(1 – 0.99)/log(1 – *p*_{s})⌉, where *p*_{s} is the corresponding success probability from Fig. 4A. The total computation time required to obtain ground states (100% accuracy) grows rapidly with problem size *N*, although it is still small in absolute terms for *N* = 100: less than 200 ms. The growth in total computation time is far less severe when the required solution accuracy is reduced; for 95% accuracy, the required time increases from ~1 ms when *N* = 10 to ~4 ms when *N* = 100; for 90% accuracy, the required time remains roughly constant at ~1 ms.

The OPO amplitudes and cut sizes/energies are shown in Fig. 4, C and D, respectively, as a function of the number of roundtrips for a single run on a single nonregular random graph instance with *N* = 100 vertices and edges. We see that the global optimum is reached within 100 μs, and the system reaches a steady state after ~120 μs. This graph instance is one of 10 instances that were generated having 495 edges. We generated 10 random graphs with a fixed number of edges, for each of , and ran these instances on our apparatus. The results of this study of the performance on random graphs as a function of edge density are shown in Fig. 4E; the system was able to find good solutions for all the attempted densities.

OPO Ising machines have the potential to harness a number of quantum features, including squeezing and superposition (*20*). Pulsed OPOs already have a substantial quantum-mechanical nature (*28*), and networks of OPOs can generate spatial multimode entanglement (*19*). The experimental system reported in this work is well-described by a quantum-mechanical model (*20*) [with very good agreement between our experimental results and numerical simulation results (*26*)]. However, the extent to which classical models can capture the behavior of OPO Ising machines is not yet known, and follow-up studies are needed to determine the fundamental computational power of OPO Ising machines. Two future modifications to the experiment that could increase the relevance of quantum mechanics are reducing the cavity round-trip loss and injecting squeezed vacuum states into the open port of the measurement beamsplitter (*29*), which is predicted to improve the success probability (*26*).

Although we find the overall direction of OPO Ising machines promising, the techniques we have demonstrated are not necessarily restricted to OPO Ising machines, and photonic-AQC and boson-sampling (*30*) experiments, among others, may be able to adapt our methods.

## Supplementary Materials

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

Materials and Methods

Supplementary Text

Figs. S1 to S5

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Materials and methods are available as supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**This research was funded by the Impulsing Paradigm Change through Disruptive Technologies (ImPACT) Program of the Council of Science, Technology and Innovation (Cabinet Office, Government of Japan). P.L.M. was partially supported by a Stanford Nano- and Quantum Science and Engineering Postdoctoral Fellowship. We thank K. Leedle, A. Ceballos, K. Wen, and Z. Wang for technical assistance and M. Digonnet, B. Lantz, T. Onodera, E. Ng, T. Leleu, C. Limouse, D. Gray, G. Tabak, and N. Tezak for helpful discussions.