Gaussian Optical Ising Machines

TL;DR

This paper demonstrates that simplified Gaussian optical Ising machines can effectively find low-energy states with success probabilities comparable to more complex systems, enabling larger and more practical implementations.

Contribution

It introduces simplified Gaussian dynamics for optical Ising machines and proposes designs without parametric gain that maintain high performance.

Findings

Gaussian dynamics suffice for high success probabilities

Simplified designs remove the need for parametric gain

Potential for larger, scalable optical Ising machines

Abstract

It has recently been shown that optical parametric oscillator (OPO) Ising machines, consisting of coupled optical pulses circulating in a cavity with parametric gain, can be used to probabilistically find low-energy states of Ising spin systems. In this work, we study optical Ising machines that operate under simplified Gaussian dynamics. We show that these dynamics are sufficient for reaching probabilities of success comparable to previous work. Based on this result, we propose modified optical Ising machines with simpler designs that do not use parametric gain yet achieve similar performance, thus suggesting a route to building much larger systems.

Equations

$$H=-\sum _{i,j}{J}_{ij}{\sigma }_i{\sigma }_j$$

$$\to \text{d}+n\eta \text{J}\text{c}$$

$$\to \text{S}\text{M}{\text{S}}^T$$

$$\to \text{S}\text{d}$$

$$\to \text{T}\text{M}{\text{T}}^T+\text{R}\text{V}{\text{R}}^T$$

$$\to \text{T}\text{d}$$

$$\left(\begin{array}{cc}\text{A}& \text{B}\\ {\text{B}}^T& \text{C}\end{array}\right)$$

$$=\text{A}-\text{B}(\text{P}\text{C}\text{P}{)}^{+}{\text{B}}^T$$

$$={\text{d}}_A+\text{B}(\text{P}\text{C}\text{P}{)}^{+}(\text{c}-{\text{d}}_C)$$

$$\to \text{C}\text{M}{\text{C}}^T+\text{S}\text{V}{\text{S}}^T$$

$$\to \text{C}\text{d}$$

Full text

Introduction

Combinatorial optimization problems, such as the travelling salesman problem, appear in many disciplines cook1998combinatorial . However, finding an optimal solution to combinatorial problems is a hard task for conventional computers. Special-purpose hardware that can solve such problems more efficiently than conventional computers is therefore an active area of research. Examples include analog electrical circuits tank1986simple ; yoshimura2015uncertain , molecular computing adleman1994molecular ; lipton1995dna , and more recently adiabatic quantum computing kadowaki1998quantum ; farhi2001quantum .

One combinatorial optimization problem of particular note is the Ising problem barahona1982computational . The Ising problem, in the absence of an external magnetic field, consists of finding the configuration of a network of coupled spins that minimizes the following Hamiltonian:

[TABLE]

where $\sigma_{i}$ and $\sigma_{j}$ are the values of the spins of sites $i$ and $j$ that can be either -1 or 1, and $J_{ij}$ is the $(i,j)$ entry of a matrix J describing the spin-spin couplings. The Ising problem maps onto several physical and combinatorial problems, such as the maximum cut problem karp1972reducibility . It is known to be NP-hard, but special-purpose hardware may find solutions to the problem faster than conventional computers.

It has recently been suggested that a train of coupled optical pulses in a cavity undergoing parametric amplification could be used as such special-purpose hardware wang2013coherent ; haribara2016coherent . In their final state, these pulses oscillate as optical parametric oscillators (OPOs) with either a [math] or a $\pi$ phase with respect to the pump light, and these two phases can be used to encode up or down spin directions. Coupling between pulses can be arranged in such a way that the system preferentially oscillates in a configuration that minimizes the Hamiltonian in equation 1.

This coupling can be achieved in an optical delay line (ODL) architecture, where all the pulses coherently interfere with each other. ODL Ising machines have been experimentally demonstrated marandi2014network ; inagaki2016large , although with low connectivity between the encoded spins. Increasing connectivity and scale is an active area of research hp2017 . More recently, a new measurement and injection feedback (MIF) scheme, that implements the required coupling using partial homodyne measurements and electronic feedback, was demonstrated with up to two thousand coupled pulses inagaki2016coherent ; mcmahon2016fully . These devices were shown to outperform several other combinatorial optimization algorithms.

In light of their experimental demonstrations of OPO Ising machines, Inagaki et al. inagaki2016coherent and McMahon et al. mcmahon2016fully , raised interesting questions concerning the nature and role of quantum features and thresholding behaviour in the operation of their devices. Given the complex dynamics of the OPO and pump fields below and above threshold, elucidating the relevant computational mechanisms presents challenges.

In this work we study MIF and ODL Ising machines that operate with simplified dynamics described by the Gaussian state formalism adesso2014continuous . We show that these simplified dynamics are sufficient to attain high success probabilities similar to what has been experimentally demonstrated. Based on these results, we propose simplifications to experimental realizations of these devices, which we envisage will be of use to scaling them up to many more coupled pulses.

Theoretical Framework

Figure 1 provides an overview of the two different types of optical Ising machines considered in this work. In both cases, a train of pulses, initially in a vacuum state, circulates in a fiber loop with parametric gain. These pulses are coupled with each other via either MIF (Fig. 1a) or ODL (Fig. 1b). In the MIF scheme, this coupling is produced by splitting off a fraction of each pulse, upon which homodyne measurement is performed. Each pulse is then displaced in phase space via optical feedback, with an amplitude and phase determined by the collective measurement results of all the other pulses. In the ODL scheme, a fraction of each pulse is repeatedly picked off, delayed, and then interfered with a consecutive pulse. In both cases, the gain or the feedback is increased with each round trip and the system settles into a final configuration of steady state OPO pulses with well-defined phases. The sign of the phases of the pulses can then be measured and mapped onto spin orientations in an Ising model, in which the spin-spin coupling is determined by the optical coupling. The system preferentially settles into a configuration that corresponds to a low energy in the Ising model.

OPO Ising machines are usually described using the stochastic master equation formalism drummond1981non ; maruo2016truncated , which provides a picture of the interaction between the pump and the signal both below and above threshold. In this work, we consider optical Ising machines that are governed by simplified dynamics described by the Gaussian state formalism adesso2014continuous . The Gaussian state formalism applies to optical states that have a Gaussian quasi-probability distribution in phase space, such as vacuum, coherent states, and squeezed states. An $N$-mode Gaussian state is fully described by its $2N\times 2N$ covariance matrix M and $2N$-length displacement vector d that characterize the Gaussian quasi-probability distribution. A Gaussian process maps a Gaussian state to another Gaussian state, and is described by simple matrix operations on the covariance matrix and displacement vector.

We note that in OPO Ising machines, many physical processes are Gaussian processes. These include loss, beam splitting, displacements, and homodyne detection. Parametric amplification can also be approximated as a Gaussian squeezing operation on the signal mode as long as pump depletion is negligible, for example in an OPO operated below threshold milburn1981production . A non-Gaussian model is necessary to account for the full interaction between the pump and the signal; this interaction is not considered in our work.

Gaussian Measurement-feedback Ising machines

We simulate a Gaussian analogue of the MIF Ising machine demonstrated by McMahon et al. mcmahon2016fully . To provide a comparison of our results to their work, we calculate the success probabilities for finding the ground state energies of specific Ising systems, where the spin-spin couplings are given by Möbius ladder graphs (see figure 2a). The $(i,j)$ entry of the J matrix corresponding to these Möbius ladders is $-1$ if nodes $i$ and $j$ are connected in the graph representing the Möbius ladder, and [math] otherwise. McMahon et al. found that their experimental system could find the correct ground state energy with a probability going from nearly 100% for small systems to about 20% for a 100-spin system.

Our simulation proceeds as follows. The pulses start in a vacuum state. In every loop, each pulse undergoes squeezing with a squeezing parameter of 0.2, followed by 30% loss. We then pick off 10% of each pulse on a beam splitter, on which we perform a simulated homodyne measurement. This is done by randomly drawing an array of numbers c from the marginal Gaussian probability distribution determined by the covariance matrix and displacement vector of the measured light (see Methods). As described by McMahon et al. we then displace the modes as follows:

[TABLE]

where $n$ is the loop number and $\eta=0.001$ is the feedback strength. With these parameters, the pump power is initially roughly 25% below the oscillation threshold of the system with $\eta=0$, which is similar to what was used in the experiment by McMahon et al. We note that when $\eta>0$, the feedback term in equation 2 brings the system above threshold.

After 300 loops through the system, we take the sign of the elements of d to be the orientation of the corresponding spins, with which we calculate the energy of the spin configuration using equation 1. We repeat this process with graphs of up to 120 spins. The probabilities of finding the correct result for these different sizes, estimated from 300 trials for each size, are shown in figure 2a. We find a probability greater than $50\%$ of finding the correct result for sizes up to about 90 spins.

The simulation above does not include gain saturation that would be present in a realistic machine, for example due to pump depletion, which could affect the success probability. To investigate the influence of saturation, we repeat our MIF simulations with the feedback displacement capped at 1000. In the case of a Möbius ladder graph the performance is adversely affected (Fig. 2A). However, we find that the change in performance is different for different graphs. Following McMahon et al. mcmahon2016fully , we calculate the success probability for 50 randomly chosen 16-spin cubic graphs, for which each spin is connected to exactly three others, using 100 trials per graph. Figure 3 compares our results for a saturated and an unsaturated MIF Ising machine. With saturation, low success probabilities are not observed. However, not every graph benefits. The success probability is only increased in 28 out of 50 cases. We observe similar outcomes for different values of the maximum feedback displacement, ranging from 100 to $10^{6}$.

Gaussian Optical Delay Line Ising machines

We also study a Gaussian analogue to ODL Ising machines. In our simulation, the pulses circulate around the loop 100 times. In every round trip, each pulse undergoes squeezing with a squeezing parameter of 0.2. Then, for every pair of pulses $(i,j)$, if $J_{ij}$ is $-1$, 10% of pulse $i$ is picked off on a beam splitter, given a $\pi$ phase shift, and interfered with pulse $j$ on a 90:10 beam splitter. After 100 round trips, we obtain a full description of the final optical state in phase space within the Gaussian state formalism. To simulate multiple experimental trials, we perform 1000 homodyne measurements on this state, from which we calculate an energy based on equation 1. We use these trials to estimate the probability of success (Fig. 2b). We find a success probability greater than $50\%$ for sizes up to about 110 spins, which is similar to what we found for the MIF scheme.

Interpretation

Our simulation of Gaussian MIF Ising machines can be understood as implementing a classical discrete random walk in optical phase space, in which the position of the walker is specified by the displacement vector. Each step taken by the walker is directed by the random outcome of the homodyne measurement via the feedback. As the feedback strength increases over time, the walker tends to move towards a specific direction away from the origin. This direction corresponds to a spin configuration which minimizes the Ising energy. This process is reminiscent of other Monte Carlo optimization schemes such as simulated annealing kirkpatrick1983optimization , in which a classical random walk probes parameter space.

In our simulation of Gaussian ODL Ising machines, we deterministically generate a probability distribution which is a multivariate Gaussian function. A specific spin configuration is found by drawing a sample from this distribution. This procedure is reminiscent of boson sampling with Gaussian states and Gaussian measurements aaronson2011computational .

We note that both types of devices considered in this work can be efficiently emulated by classical computation. This can be seen as a consequence of the fact that Gaussian optical Ising machines involve only Gaussian measurements on Gaussian states bartlett2002efficient .

Simplified Ising Machines

The OPO Ising machines demonstrated by Inagaki et al and McMahon et al can potentially scale up to tens of thousands of coupled pulses haribara2016coherent . Here we propose simplified Ising machines that may provide a route to building large scale devices.

We first consider a Gaussian MIF Ising machine in which we simply remove the squeezing but maintain the measurement and feedback. The pulses in the system are therefore coherent states. In the absence of the gain that previously partially compensated for the loss, we double the feedback strength to $\eta=0.002$. The random walk process implemented by these simulations follows the same pattern as in the simulation with parametric gain, but random numbers are now drawn from a narrower distribution around the average displacement due to the absence of squeezing. Even though the phase of a coherent state is not intrinsically bistable as for OPOs, the use of measurement and injection that displaces along only one quadrature axis effectively makes the phase bistable. Our simulation results are shown in figure 2a. We find that removing the squeezing has only minor influence on the success probability. We also find, as in the case with squeezing, that the success probabilities for different graphs can be modified by saturated feedback.

Another alternative to an OPO-based MIF Ising machine involves the use of phase insensitive gain media instead of parametric gain, as has been proposed utsunomiya2011mapping . Such a scheme would allow for tuning the uncertainty on the quadratures of the pulses, thus adding an additional controllable degree of freedom to the random walk implemented by MIF Ising machines.

We also simulate the Gaussian ODL Ising machine as before, except the squeezing operation is replaced with phase insensitive gain. We note that in this case, the phases of the pulses are no longer bistable, but they remain correlated due to their mutual couplings. Our results are shown in figure 2B. The performance with phase insensitive gain is essentially identical to that with squeezing.

Conclusion

We introduce Gaussian optical Ising machines, and find that Gaussian dynamics are sufficient for finding low-energy solutions of the Ising problem, with performance comparable to demonstrated OPO based machines. Based on this finding, we propose simplified Gaussian approaches that suggest a route to developing future large scale optical Ising machines.

We thank Jan Sperling, Alireza Marandi, Peter McMahon, Edwin Ng and Tatsuhiro Onodera for helpful discussions. We acknowledge support from NWO Rubicon, the European Research Council, the UK Engineering and Physical Sciences Research Council (project EP/K034480/1 and the Networked Quantum Information Technology Hub), from The John Templeton Foundation, and from the European Commission (H2020-FETPROACT-2014 grant QUCHIP).

Methods

The following provides more detail as to how we implemented our simulations using the Gaussian state formalism. In our simulations, we used the convention that the covariance matrix of vacuum is $\frac{1}{2}\mathbb{1}$, but the following formulae do not depend on the choice of convention.

.1 Simulation of Squeezing and Loss

Single mode squeezing on mode $i$ of a $N$-mode Gaussian state changes its covariance matrix M and displacement vector d as follows:

[TABLE]

where S is unity except for entries $2i-1$ and $2i$ along the diagonal which are equal to $e^{r}$, where $r$ is the (real) squeezing parameter.

Optical loss on mode $i$ is modelled by a beam splitter with transmission $t$ acting between mode $i$ and an ancilla vacuum mode. This changes M and d as follows:

[TABLE]

where T is unity except for entries $2i-1$ and $2i$ along the diagonal which are equal to $\sqrt{t}$, R is unity except for entries $2i-1$ and $2i$ along the diagonal which are equal to $\sqrt{1-t}$, and V is the covariance matrix for an $N$-mode vacuum.

.2 Simulation of Homodyne Measurements

A homodyne measurement on the $N$-th optical mode of an $N$-mode Gaussian state described by covariance matrix M and displacement vector d is simulated as follows. M can be written in block-diagonal form:

[TABLE]

where C is the reduced $2\times 2$ covariance matrix for the homodyned mode, A is the reduced covariance matrix for the remaining $N-1$ modes, and B keeps track of the correlations between these modes. The homodyne measurement projects mode $N$ onto an infinitely squeezed state with a displacement vector c, which is drawn from the projection of the marginal Gaussian probability distribution described by C and d onto the chosen measurement axis. The remaining $N-1$ modes are then described by a new covariance matrix A’ and displacement vector d’:

[TABLE]

where P is a projector onto the measured quadrature of mode $N$, + denotes the Moore-Penrose pseudo-inverse, and $\textbf{d}_{A}$ and $\textbf{d}_{C}$ are the original reduced displacement vectors for the first $N-1$ modes and for the $N$-th mode, respectively.

.3 Simulation of Phase Insensitive Gain

Phase-insensitive gain on a given mode $i$ of an $N$-mode state can be modelled in the Gaussian state formalism by using an ancilla vacuum mode, implementing a two-mode squeezing operation with squeezing parameter $r$ on these two modes, and tracing out the ancilla mode. In the $(x_{1},p_{1},...,x_{N},p_{N})$ basis in optical phase space, the covariance matrix M and the displacement vector d for the amplified mode are thus transformed as follows:

[TABLE]

where C is unity except for entries $2i-1$ and $2i$ along the diagonal which are equal to $\cosh(r)$, S is unity except for entries $2i-1$ and $2i$ along the diagonal which are equal to $\sinh(r)$, and V is the covariance matrix for an $N$-mode vacuum. In our simulation of an ODL Ising machine with phase-insensitive gain, we used a squeezing parameter of $r=0.6$, which yields an amount of gain $\cosh(r)\approx 1.19$ that is similar to the gain $e^{r_{1}}\approx 1.22$ provided by single mode squeezing with $r_{1}=0.2$.