A fully programmable 100-spin coherent Ising machine with all-to-all connections

See allHide authors and affiliations

Science  04 Nov 2016:
Vol. 354, Issue 6312, pp. 614-617
DOI: 10.1126/science.aah5178

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.

Science, this issue pp. 603 and 614


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 (35). 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 (1114), 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 (1720), 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 Embedded Image (i = 1, ..., N) that minimizes the energy function Embedded Image, where the particular problem instance being solved is specified by the N × N matrix J (with elements Jij) and the length-N vector h (with elements hi).

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

Fig. 1 Experimental schematic of a measurement-feedback–based coherent Ising machine.

A time-division–multiplexed pulsed degenerate optical parametric oscillator is formed by a nonlinear crystal [periodically poled lithium niobate (PPLN)] in a fiber ring cavity containing 160 pulses. A fraction of each pulse is measured and used to compute a feedback signal that effectively couples the otherwise-independent pulses in the cavity. IM, intensity modulator; PM, phase modulator; LO, local oscillator; SHG, second-harmonic generation; FPGA, field-programmable gate array.

We have focused our experiments on Ising problems on undirected, unweighted graphs (V, vertices; E, edges), where Jij = –1 when spin i and spin j are connected, and Jij = 0 otherwise, and for which the linear (Zeeman) terms hi are zero. An Ising problem of this form is equivalent to the Max-Cut problem on the underlying graph (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,...,N for an Ising problem is given by Embedded Image. The same spin configuration {σi} represents a cut of size Embedded Image; 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.

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 (ci) 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 Tcomp = TrtNrt, where Trt = 1.6 μs is the single round trip time and Nrt is the number of round trips. Each spin σi in the Ising problem is represented by a single OPO pulse; if the in-phase component ci of the ith OPO pulse is less than zero (ci < 0), we make the spin assignment σi = –1, and if the in-phase component is greater than zero (ci > 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 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 Nrt = 300 round trips, a computation time of Tcomp = 480 μs.

Fig. 2 Results with N = 16 cubic graphs.

(A) A Möbius ladder graph with N = 16 vertices. (B) The evolution of the in-phase components ci of the N = 16 OPO pulses as a function of the computation time. (C) The graph cut size achieved as a function of the computation time. (D to F) Histograms of obtained solutions in 100 runs for the graphs shown in the insets. (G) Histogram of the observed probabilities of obtaining a ground state in a single run, for all 4060 unweighted, undirected cubic graphs with N = 16 vertices.

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

Fig. 3 Results with various-size Möbius ladder graphs.

(A) Observed probability of obtaining a ground state of the Möbius ladder graph in a single run, as a function of the size N of the graph. Multiple 100-run batches were performed for each graph size to obtain the standard deviations, which are shown as error bars. (B to D) Histograms of obtained solutions in 100 runs for the graphs shown in the insets.

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.

Fig. 4 Results with various-size and various-density random graphs.

(A) Observed probability of obtaining a solution whose cut size is at least x% of the global optimum (maximum cut), as a function of graph size N, for random cubic graph instances. Error bars indicate 1 SD, which is dominated by the difference in difficulty between the various problem instances. (B) The runtime that would be required to obtain a solution of a particular accuracy with 99% probability. (C) The evolution of the in-phase components ci of the N = 100 OPO pulses as a function of the computation time, for a single run with the graph shown in the (D) inset. (D) The graph cut size achieved as a function of the computation time. (Inset) The graph being solved. (E) Observed success probability of obtaining a solution with a particular accuracy as a function of the density of edges in the graph. Experiments were performed on randomly generated N = 100-vertex graphs with fixed numbers of edges. Error bars indicate 1 SD.

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 Tcomp × ⌈log(1 – 0.99)/log(1 – ps)⌉, where ps 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 Embedded Image 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 Embedded Image, and ran these instances on our apparatus. The results of this study of the performance on random graphs as a function of edge density Embedded Image 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

Materials and Methods

Supplementary Text

Figs. S1 to S5

References (3158)

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. 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.

Stay Connected to Science

Navigate This Article