System and method for efficient sampling of ground states and low energy Ising spin configurations using a coherent Ising machine

Abstract

A system and method for efficient sampling of ground states and low energy Ising configurations. The system may be implemented using the nonlinear stochastic dynamics of a Measurement Feedback Based Coherent Ising Machine (MFB-CIM). The discrete-time Gaussian state model of the MFB-CIM can capture the nonlinear dynamics. The system and method require many fewer loud trips to sample than in other known systems.

Description

[Technical field]

[0001] This application claims priority from U.S. Provisional Application No. 63 / 157,673, filed March 6, 2021, the entirety of which is incorporated by reference.

[0002] Appendices A and B (2 pages) have details on the quadratic equations or motions for crystal propagation and on evaluating expectation values ​​of the quadratic operators.

[0003] Each of the above Appendices forms part of and is incorporated herein.

[0004] The present disclosure relates generally to systems and methods for efficient sampling of ground-state and low-energy Ising spin configurations, and more particularly to systems and methods that use coherent Ising machines to implement efficient sampling of ground-state and low-energy Ising spin configurations. [Background technology]

[0005] For decades, the Ising model has served as an important conceptual bridge between the fields of physics and computation. A host of important combinatorial optimization problems have efficient mappings to the problem of finding ground states of the Ising model [1], while the simple and very general form of the model implies that Ising-like interactions are ubiquitous across a diverse array of systems [2]. Formally, the Ising model is defined by the Ising Hamiltonian:Σ i≠j J ij σ i σ j A set of spins σ with configuration energy given by i=±1, and in general the Ising problem of finding the spin configuration that minimizes this energy is currently intractable for classical computers [3]. As a result, there has been considerable interest in leveraging physical Ising-like systems as dedicated computational hardware to tackle problems such as combinatorial optimization, with ongoing research on platforms ranging from quantum annealers built from microwave superconducting circuits [4, 5] to coherent Ising machines based on networks of nonlinear optical oscillators [6, 9], among many others [10, 15].

[0006] However, while combinatorial optimization is often focused on finding just one of the ground-state Ising spin configurations, it is desirable in many applications to obtain many or all degenerate ground-state configurations, and in some cases to also sample many low-energy configurations

[16] . Such sampling capabilities are particularly useful for applications that involve obtaining distributed information about the spin configurations within an Ising model, such as estimating the ground-state entropy of physical simulations with Ising-like interactions or implementing Boltzmann machines as generative models for machine learning [17, 19]. In industrial settings, accessing a pool of candidate solutions to an optimization problem can make the process more efficient and flexible; for example, in drug discovery [20, 23], structure-based lead optimization can generate a large number of candidate molecules for simultaneous testing. Recently, it has also been pointed out that when decomposing a large-scale optimization problem into subproblems that are solved separately (e.g., to address hardware limitations), a better solution to the original problem can be constructed using multiple low-energy samples, rather than just being optimal for each subproblem [24, 25]. However, Ising solvers designed for combinatorial optimization are not necessarily suitable for sampling all ground states and / or low-energy states. For example, commercial quantum annealers by D-Wave Systems have shown success in finding ground states of Ising problems, and their operating principles can result in an exponential bias in the distribution of degenerate ground state samples [26, 28]. It is therefore desirable, and it is for this reason that the present disclosure is directed, to provide a system and method for efficient sampling of ground states and low-energy CIM states that overcomes the above limitations and problems of conventional techniques. [Brief description of the drawings]

[0007] [Figure 1]FIG. 1 illustrates a measurement-feedback-based coherent Ising machine (MFB-CIM). [Figure 2A] FIG. 1 illustrates how quantum noise can affect the nonlinear dynamics of CIM. [Figure 2B] FIG. 1 illustrates how quantum noise can affect the nonlinear dynamics of CIM. [Diagram 3] FIG. 13 shows the trajectories of the discrete-time model for different cavity decay times. [Figure 4] FIG. 1 illustrates MFB-CIM sampling for low-energy spin configurations for a particular N=15 Ising problem. [Diagram 5] FIG. 5 illustrates the variation of the maximum required sampling time for the spin configuration shown in FIG. [Figure 6] FIG. 14 shows the scaling of the time required to sample a single Ising ground state as a function of cavity damping. [Figure 7] FIG. 13 illustrates the required sampling time for an alternative MFB-CIM sampling model. [Figure 8] FIG. 13 illustrates the scaling of the sampling performance of the negative pump MFB-CIM with respect to problem size N. [Figure 9] FIG. 13 illustrates sampling performance using various alternative CIM sampling methods. [Figure 10A] 1 is a flow chart of an iterative method for efficient sampling. [Figure 10B] 1 is a flow chart of an iterative method for efficient sampling. [Figure 11] FIG. 2 illustrates a computer system that can incorporate the CIM of FIG. 1 and perform efficient sampling. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0008] The present disclosure is particularly applicable to systems and methods for efficiently sampling ground states and low energy Ising spin configurations using the disclosed measurement feedback based coherent Ising machines ("MFB-CIMs"), and it is in this context that the present disclosure will be described. However, it will be recognized that the systems and methods have greater utility as they may be implemented using alternative embodiments and implementations, which are discussed in more detail below.

[0009] Systems and methods for efficient sampling of ground states and low-energy spin configurations of Ising machines can be implemented using a specifically constructed measurement feedback based coherent Ising machine (MFB-CIM) [7, 9, 29]. The MFB-CIM is a hardware platform originally conceived for performing Ising optimization using a network of degenerate optical parametric oscillators (DOPOs) subjected to a real-time measurement feedback protocol that encodes Ising interactions within the network dynamics. In particular, the method uses a Gaussian state model to consider how quantum noise arises within the dynamics of the MFB-CIM and to address whether such stochastic nonlinear dynamics can facilitate efficient sampling of low-energy Ising configurations. While fully quantum processing of MFB-CIMs is possible

[29] , numerical investigations of large-scale systems suitable for combinatorial optimization / sampling are only possible up to the Gaussian state regime, and quantum correlations are only considered up to second order (i.e., up to the covariance of observables)

[30] . This Gaussian state approximation is consistent with all experimental MFB-CIM operating regimes known to date, while still providing an accurate treatment of important quantum-noise-driven phenomena such as squeezing / antisqueezing and measurement uncertainties and backaction [29, 31, 32], which are crucial for investigating the sampling performance but are typically neglected in mean-field models.

[0010] The potential of MFB-CIM to efficiently generate samples of degenerate basis and low-energy excited spin configurations was recently pointed out in

[33] using a Gaussian-state quantum model

[34] formulated in continuous time. In the system and method, the continuous-time model correctly captures the dynamics of MFB-CIM in the high-finesse limit, where the cavity decay time of its constituent DOPOs dominates all other system timescales. On the other hand, the inherently high bandwidth of low-finesse cavities can, at least in principle, be leveraged to significantly reduce the computational runtime; in practice, most experimental implementations of CIM (both optically coupled and measurement feedback based) utilize DOPOs operating in the low-finesse regime with short cavity decay times [6, 9]. Low-finesse systems are more conveniently described in discrete time, where the dynamics occurs by a sequence of discrete operations with respect to the system states. Theoretically, such a discrete-time model of MFB-CIM has been previously investigated in references [35, 36]. While the latter investigation used a non-Gaussian model for the quantum state, which is only numerically tractable for small problem sizes, the former effort used a Gaussian state model to investigate the linear dynamics of MFB-CIM. However, in this Gaussian model, nonlinear gain saturation, which in principle can play a non-negligible dynamic role in MFB-CIM near and above threshold, was only considered phenomenologically. To circumvent these limitations of known techniques, the system and method uses a discrete-time Gaussian state quantum model featuring a physical model for nonlinear gain saturation, which can investigate low and intermediate finesse MFB-CIM below, through and above threshold. Gaussian Quantum Model for Measurement Feedback Based CIM (MFB-CIM)

[0011] Conceptually, a coherent Ising machine (CIM) is a system of N degenerate optical parametric oscillators (DOPOs), which are nonlinear optical oscillators that exhibit saturable phase-sensitive gain. When pumped below its oscillation threshold, the state of the DOPO is well described by a quadrature-squeezed state, whereas further above the threshold, nonlinear saturation of the gain due to pump depletion stabilizes the system to be in one of two phase bistable bright coherent states (called 0 and π phase states). To encode the Ising spin within the DOPO network, the method associates these bistable phase states with the Ising spin = ±1. By engineering the interactions between the DOPOs, the system has system dynamics governed by the desired Ising coupling matrix J.

[0012] FIG. 1 illustrates at a schematic level a measurement feedback based coherent Ising machine (MFB-CIM) that may be used in the disclosed system. In this MFB-CIM, degenerate optical parametric oscillators (DOPOs) are realized as synchronously pumped and time multiplexed "signal" pulses in a single optical cavity, with the pulse interactions mediated by a synchronous real-time measurement feedback protocol. The device consists of a system of DOPOs and a feedback control module, as shown in FIG. 1. The OPOs are coupled to a second order (χ (2) ) is pumped through a medium with nonlinearity using a second harmonic pulse of light. Homodyne measurement of the intracavity pulse produces a state ρ (h) This is done by extracting the intracavity pulse by a beam splitter action with an incoming pulse having a . The measured homodyne signal is fed to an FPGA which performs a matrix-vector product. Finally, this output is used in the OPO to apply a displacement to the intracavity pulse. The Gaussian quantum model for the MFB-CIM model models the physics of this device in discrete time by providing a theoretical description of all the physical processes, labeled with Roman numerals, in the CIM.

[0013] In the measurement feedback based CIM (MFB-CIM) shown in Figure 1, the signal pulses are spaced at time intervals of 1 / f rep Therefore, to fit N pulses, the cavity length is ≈Cn / f rep For each round trip through the cavity, the signal pulse is coupled in parallel with a synchronized externally injected pump pulse to produce a nonlinear χ (2) The signal pulses are then successively tapped off through an output coupler. The output is then measured on a q-quadrature homodyne detector, providing an indirect and weak measure of the q-quadrature amplitude of the internal signal pulse. Importantly, the sign configuration of the N homodyne measurements, e.g. (sgn ω 1 , …, sgn ω N ) is the corresponding σ i←→ sgn ω i Finally, to perform the pulse-pulse interaction, the FPGA in Fig. 1 receives the homodyne result and generates an analog feedback signal υ i ∝Σ j J ij ω j and the signal is converted (for example, using a synchronous optical modulator) into i is applied to the corresponding i-th signal pulse by synchronous external injection of a feedback pulse with magnitude and phase determined by: The interference between the injected pulse and the internal signal pulse steers the system towards a low-energy Ising spin configuration, thus dynamically realizing the Ising coupling matrix J in the MFB-CIM system.

[0014] The consequence of embedding the structure of Ising couplings into the system dynamics is that the evolution of the state is governed by the interplay of nonlinearities that drive the signal amplitude towards bistable spin values; linear couplings that drive the system towards a collective configuration that minimizes the Ising energy; and quantum noise that arises from the inherent uncertainty of weak homodyne measurements and drives the system by both measurement backaction and feedback injection. As conceptually shown in Figure 2, the dynamic evolution of the signal amplitude is stochastic and nonlinear, with the preferred code configuration governed by the Ising coupling matrix. In the continuous-time limit, a convenient and intuitive picture is to think of such stochastic trajectories as noisy gradient descent on a potential landscape: as the state evolves stochastically in time, the instantaneous potential seen by any given spin also changes dynamically, guiding the system around and across local minima that arise from the interplay between nonlinearities and couplings.

[0015] An open dissipative bosonic particle system such as the MFB-CIM, which has a relatively weak single-photon induced nonlinearity and undergoes continuous homodyne measurements, can usually be well approximated by a Gaussian state. Formally, it has a Wigner function whose quantum state is approximately Gaussian, and as a consequence can be fully characterized by simply specifying a set of mean-field amplitudes and a set of covariances that describe the quantum correlations and their amplitude uncertainties. Another very useful simplification applied to the MFB-CIM is that the pulses, while interacting through measurement feedback, are still in a non-entangled state, resulting in zero covariance between different signal pulses, since the physics in Fig. 1 involves local operations and classical communication (LOCC) between the signal pulses.

[0016] Many of the operations involved in FIG. 1, including out-coupling, homodyne measurements, and feedback injection, are linear operations and have a particularly simple description when the method uses a discrete map formulation for the dynamics in which, for Gaussian states, the signal pulse is assumed to undergo a discrete input-output transformation as it passes through each optical component. This discrete map approach contrasts with the more traditional continuous-time model in quantum optics, where the amplitudes evolve continuously in time under a set of Liouvillian superoperators that describe the action of the effective system Hamiltonian as well as the losses and measurements. Of course, the time-multiplexed and pulsed nature of the setup described in FIG. 1 naturally lends itself to a discrete map model when the pulse widths are short compared to their separation. Nevertheless, the continuous-time model can be a good approximation when the cavity finesse is high: in this limit, the single round-trip gain, loss, and measurement strength are all small, resulting in small changes in the cavity state for every round-trip, so that the overall system dynamics is well described by a continuous-time differential equation. To investigate the effect of varying the cavity finesse on the sampling performance, the method however uses the more general framework of a discrete-time formulation in which the linear loss can be specified independently of other parameters.

[0017] Importantly, one operation in MFB-CIM that does not easily lend itself to a discrete-time model is the propagation of a signal pulse through the crystal. Below threshold, this operation can be well described by a linear quantum-limited gain along the q-orthogonal

[31] , which is modeled as a discrete-time squeezing operation in ref.

[35] . On the other hand, when the DOPO is near or above threshold, the forward pump pulse is depleted, which saturates the gain and leads to nonlinear dynamics. The main contribution of the following model is to prescribe an efficient numerical treatment of this latter gain saturation physics, allowing the discrete-time model to be extended into the above-threshold regime while remaining consistent with the standard continuous-time quantum optical model for MFB-CIM in the high-finesse limit.

[0018] A formal description of the Gaussian state model is given below, culminating in a fully iterative algorithm for simulating the MFB-CIM in discrete time. This disclosure also shows that in the high finesse limit, the dynamics of this model reduces exactly to the dynamics of the standard continuous-time quantum model used to process the MFB-CIM in the Gaussian state regime. Discrete-time Gaussian quantum model

[0019] In this disclosure, the time multiplexed MFB-CIM is

number

number

number

number

number

[0020] When a system is in a Gaussian state, it is completely determined only by the mean vector and the covariance matrix Σ; i.e., the quantum state is

number

number

number

number

number

number

number

number

number

number

[0021] More generally, the state

number

number

number

number

number

number

[0022] With the above basic formulation established, there are several linear operations that are necessary for the operation of the MFB-CIM before moving to the nonlinear map. These basic operations consist of a beam splitter to model losses and out-coupling, coherent injection to model feedback, and homodyne measurements. The field exchange amplitude r (i.e., the power exchange ratio r 2 ) two-mode state

number

number

number

number

[0023] Displacement into mode a

number

number

number

number

number

number

number

number

number

number

number

number

number

number

[0024] These linear maps describing the light extraction, measurement, and feedback injection are fairly simple, but dissipative linear losses also need to be addressed.

[0025] The experimental sources of loss in physical CIM vary with implementation details, but some prominent sources include crystal facet losses (due to mode-matching inefficiencies or Fresnel-reflection losses) and cavity propagation losses (due to scattering from mirrors or mode-matching inefficiencies while coupling to / from the fiber). Because crystal facet losses generally dominate in realistic experimental implementations, the method assumes for simplicity that all loss mechanisms can be lumped together and applied by a pair of partial beam splitters placed before and after the crystal. Like the out-couplers used for measurements, these beam splitters extract the intracavity light, but instead of the outgoing pulse being measured by homodyning (which would cause backaction on the state), this external pulse cannot be measured, instead it simply traces it out partially, resulting in dissipation on the state. Nonlinear Crystal Propagation

[0026] The most challenging part of the discrete-time model concerns the propagation of the pulse through the nonlinear crystal, which is a dynamic non-Gaussian process, in contrast to other operations, including measurements and feedback, which can all ideally be treated as Gaussian operations. Each incoming signal pulse

number

number

number

[0027] In general, the Hamiltonian (11) can give rise to both entanglement and non-Gaussianity in the coupled state between the pump and the signal pulse, and requires a fully coupled Hilbert space of two modes to be adequately described. To fit the crystal propagation into a Gaussian formulation, we derive equations of motion (EOM) for the Gaussian moments of the pump and signal pulses generated by (11), while assuming that the non-Gaussianity of the state (characterized by higher-order moments) remains negligible. This approximation is valid when the DOPO has a large saturation photon number, i.e., a single photon only induces a small gain saturation. The EOM can then be numerically integrated from the input facet to the output facet of the crystal,

number

[0028] One simple way to restrict the quantum dynamics to a Gaussian subspace is to take the Heisenberg equations of motion generated by (11) for the orthogonal operators and perform a moment expansion up to second order

[40] .

number

number

number

number

number

number

number

number

[0029] The evolution of the first moment can be obtained simply by taking the expectation with respect to the above equation. To resolve the product, the relationship

number

number

number

number

number

number

[0030] The complete equations of motion derived under this procedure are provided in Appendix A. In general,

number

number

number

number

number

number

[0031] The final Gaussian state EOM can be numerically integrated to perform the crystal propagation map of motion on the crystal propagation map of (12). The mean-field equation is

number

number

number

number

number

number

[0032] Using the above components and transformations to model the MFB-CIM, a method 1000 (shown in FIGS. 10A-10B) for an iterative procedure to generate the dynamics of the MFB-CIM is now described.

number

number

[0033] 10A and 10B show how the iterative method shown in FIG. 10A and FIG. 10B can be used to determine the state of the i-th signal pulse.

number

number

number

number

[0034] The method can then perform / determine 1004 the crystal propagation. According to (12), the crystal propagation is described by a Gaussian map, the joint correlation signal pump states and subsequent partial traces of the pump modes:

number

number

number

[0035] The method can then determine 1006 the output facet losses. This process is exactly the same as for the input facet losses. Assuming that we lump the total system losses symmetrically between the input and output losses around the crystal, the process of

[18] can again be applied for this.

[0036] The method then performs photo-extraction and homodyne measurement (1008). The homodyne measurement consists of two sub-processes. First, a portion of the internal signal pulse is photo-extracted, which is then converted into a map.

number

number

number

[0037] The method can then inject (1010) measurement feedback as shown in FIG. 10B. To do so, the method applies a coherent displacement to the signal pulse based on the feedback signal calculated by the FPGA in the CIM to perform the Ising combination. The feedback term is:

number

number

[0038] 10B, the method determines 1012 whether there are more pulses 1012 and loops back to process 102 so that the process in the method is applied to each pulse. After the process is applied to each pulse i=1,...,N, it completes one round trip through the CIM cavity. Note that the exact sequence of the following operations will technically depend on the details of how the fiber cavity is laid out and the relative flight times between the optical components. Nevertheless, general features such as steady-state behavior should be robust to the exact choice of ordering, and if precise transient behavior is desired, the above procedure can be rearranged to more accurately model a particular cavity layout. Once the iterations are completed for all of the pulses, the method outputs the results (1014). Reduction to a continuous-time Gaussian model

[0039] A conventional continuous-time quantum model for MFB-CIM in the Gaussian regime can be compared to the above method, showing that the dynamics produced by the disclosed discrete-time model match those of the continuous-time model in the high-finesse limit. Continuous-time Gaussian quantum model

[0040] The standard approach to modeling MFB-CIM is based on input-output theory [32, 41], which describes an open quantum system weakly coupled to a set of external reservoirs. In this formulation, the dynamics is specified by a system Hamiltonian that captures the unitary evolution, and a set of Lindblad operators that describe the interaction of the system with the reservoirs.

[0041] For MFB-CIM, a system of N OPOs is

number

number

number

number

[0042] The Hamiltonian consists of two coherent effects. The first is of the form

number

number

number

number

[0043] Measurement record m j Since (t) constitutes a continuous weak measurement of the system state, the dynamics of the system is j In standard input-output theory, such dynamics is expressed by a stochastic master equation (SME)

[42] .

number

number

number

[0044] From the SME, the conditional evolution of any desired observable can be obtained. To establish a correspondence with the discrete-time model, we use the in-phase, quadrature

number

number

number

number

number

number

number

number

number

number

number

number

[0045] The continuous-time dynamics of the form (27) can be obtained from a discrete-time model in the high-finesse limit. In the high-finesse limit, each discrete operation

number

number

[0046] The parameter δ is introduced such that δ → 0 formally defines the high-finesse limit. The MFB-MIC round-trip time (measured by the wall clock) is Δt = N / f rep . The model parameters appearing in method 1000 scale as follows:

number

number

[0047] Each of the operations in method 1000 discussed above can be extended to their leading order correction in δ. Since the q- and p-orthogonal parts of the dynamics are separated, as explained in the subsection above for the discrete-time model, and as shown in (27) for the continuous-time model, we use the following

number

[0048] First, the linear loss in the facets can be given by (18). Using (5), this gives the mapping

number

number

number

number

number

number

number

number

number

number

[0049] Third, for the measurement process and the light extraction step, the signal state is:

number

number

number

number

number

number

[0050] Finally, taking into account the injection feedback according to (21) and the measurement results (35), the applied displacement is

number

number

number

[0051] All maps in a single round trip by combining effects from (29), (31), (32), (35), and (36) up to first order in δ (again, ignoring order). A discrete-time finite difference map describing one round trip limit for the continuous-time differential equation, with updated means and variances denoted by primes.

number

number

[0052] Finally, as an alternative to the above approach, in which the correspondence is derived by a Picard iteration for a (c-number) Gaussian state EOM, it is also possible to reach the same conclusions by a quantum input-output analysis of the crystal Hamiltonian (11). For example, for one propagation, the crystal can be treated as a unitary operation

number

number

number

number

[44] the above unitary over one round-trip time by writing

number

number

number

number

number

number

[0053] In the continuous-time theory over many round trips, (42b) is precisely the gain / squeezing part of the continuous-time system Hamiltonian (23), while (42c) is the continuous-time Lindblad operator in the continuous-time model.

number

[0054] Numerical simulations of the discrete-time Gaussian model are provided along with some representative trajectories of the model dynamics to define an appropriate metric for the sampling performance. Model parameters

[0055] For numerical results, it is useful to define a new set of parameters that scales the model more favorably by keeping certain qualitative features constant. The dynamics of an uncoupled classical DOPO can be calculated using four parameters: (1) the cavity photon 1 / e in the absence of pumping 2- Decay time; pump field β and its threshold β th The pump parameter r=β / β th ; and (3) the saturation photon number n sat Each of these quantities can be expressed in terms of the model parameters given above. For convenience, the method is

number

[0056] First, in the absence of pumping or nonlinearity, the number of photons decreases by 1 / e due to linear losses and light extraction. 2 The number of round trips required to attenuate by a factor of T decay teeth,

number

number

[0057] Classically, the threshold pump field is defined to be the value of β (i.e., the input pump pulse amplitude) such that the exponential gain experienced by a small signal (i.e., a signal pulse with an extinction amplitude) input into the crystal is exactly balanced with the attenuation due to linear losses and light extraction. The pump parameters are then simply determined by adjusting the pump field to this threshold β th Therefore, the regulation is

number

number

[0058] Classically, the saturation photon number is the square of the mean-field signal amplitude at steady state for r=2, with contributions from linear loss / light extraction, parametric gain, and nonlinear gain saturation. Because this feature involves nonlinear terms in the crystal EOM at finite signal amplitudes, the exact value of the saturation photon number can depend on the cavity layout for low finesse cavities. However, in the high finesse limit, its value is

number

number

number

[0059] Finally, for the injected feedback coupling with r<1, the feedback gain J required for the system to exceed the threshold with the feedback gain is 0 teeth,

number

number

[0060] Figure 3 shows some representative dynamics of MFB-CIM run on a small N = 8 problem instance. These trajectories are shown to be consistent with T decay Note especially how the θ changes as a function of θ. In all cases, θ is 10T. decay By running a simulation on the low finesse system (T decay =4) to the intermediate finesse system (T decay =16) with high finesse system (T decay = 64), the internal Gaussian state dynamics eventually converges to a single continuous-time trajectory in the high-finesse limit as shown in the rightmost column. As expected, the homodyne record w i (Thus, the measured Ising energy -Σ j J ij w i w j ) is the light extraction ratio

number

[0061] As can be seen in Fig. 3, the dynamics of the MFB-CIM drive the system towards above-threshold steady states where the sign configuration of the DOPO encodes the low-energy spin configuration of the Ising problem. At the same time, the exact spin configurations found by the MFB-CIM in these dynamics are probabilistic, driven by the measurement backaction and feedback through homodyne detection. In certain operating regimes, the MFB-CIM can be used to probabilistically sample different spin configurations, simply by launching the same system under different realizations of the measurement noise. Each "launch" of the MFB-CIM takes

number

[0062] The bar graph in Fig. 4(a) illustrates this procedure for a particular N=16 problem instance by recording for each low-energy Ising spin configuration the number of trajectories in which that configuration appeared at least once. Over 1000 trajectories, the method and system readily obtains multiple samples of all spin configurations of interest that exhibit relatively fair sampling, at least for this instance. Nevertheless, there are some systematic biases in the sampling; i.e., spin configurations at a given energy level are not necessarily sampled uniformly. These biases are generally problem-dependent, but also depend on the model parameters chosen for the sampling process.

[0063] To fully quantify the efficiency and fairness of the sampling, however, it is not sufficient to simply count the trajectories in which each spin configuration appears, since certain configurations may systematically appear later than others within any given trajectory. These differences in sampling times within trajectories are illustrated in Fig. 4(b), which shows, for each spin configuration considered in Fig. 4(a) (up to the full spin flip), a (vertical) histogram of the first time that the configuration appeared in each trajectory, if it appeared at all. In all cases, there is an initial transient period (~T decay round trip), with most of the distribution peaking within the short decay time of the transient. However, the exact distribution of this first sampling time varies from spin configuration to spin configuration, with some configurations featuring a sharper peak and others having a longer tail. As a result, spin configurations that appear less frequently per trajectory may nevertheless be efficient to sample if they appear early within the trajectory in which they do appear.

[0064] A "required sampling time" metric can be defined to take these effects into account, including bias in the overall count, transient time costs, and variations in the first sampling time distribution. For example, the homodyne record

number

number

number

number

number

number

[0065] This defined empirical measure of the required sampling time is affected by the model parameters, in particular the feedback gain and the pump strength. Figure 5 shows the maximum required sampling time T samp We show how �� varies across different model parameters for the problem instance and spin configuration considered in Fig. 4. Efficient sampling in MFB-CIM is relatively robust across a wide range of system parameters. The most important parameters are the feedback gain α and the pump parameter r, which show a sharp threshold in the sampling performance near the estimated linear threshold of MFB-CIM. The required sampling time is a function of the faster cavity decay time (i.e., T decay = 4), which is a reflection of the fact that the low-finesse cavity spends fewer round trips in the transient period and can result in a low-energy sample more quickly. There also appears to be a slight advantage to using a system with a lower escape efficiency (higher unmeasured losses), which can be attributed to the effective reduction of measurement backaction in such a regime.

[0066] One particularly interesting aspect of the sampling behavior in the MFB-CIM is that the required sampling time is T decay The goal of this study is to scale the problem with the finesse of the system, measured by . To check whether this scaling is robust with respect to the choice of problem instance, a set of integer-valued Sherrington-Kirkpatrick Ising problems is considered with range 1 (SK1), which is equivalent to the set of MAX-CUT problems with signed binary edge weights. Figure 6 shows the scaling of T over 50 SK1 problems. sampWe denote the distribution of T (σ), which for concreteness is chosen to be the first lexicographic configuration that gives the ground energy of the problem. As shown, there is a spread in the required sampling times among the problem instances, and the distribution can be nicely characterized by a mean and median that show a clear monotonic decrease with decreasing decay time. Interestingly, this scaling is consistent with the T decay It persists down to very low decay times, on the order of 1. In practice, for this problem size, performance is better than round-trip decay of T decay begins to approach 1 exponentially as a function of (e.g., R loss ~1-e -Tdecay ) point, T decay It only saturates and declines at ∼0.2. At this point, the sensitivity of the system to the system parameters precludes any further significant gains when reducing the sampling time. Continuing to improve sampling performance into the low-finesse regime is an important motivation for our development of a discrete-time Gaussian model. Sampling of alternative models of MFB-CIM

[0067] Although there has thus been an increased focus on developing a general model for MFB-CIM in the Gaussian state approximation and on investigating the dynamic role of quantum noise in its conventional operation (with parametric gain, homodyne measurement / feedback, measurement backaction, and gain saturation), it is also useful to consider alternative models or modes of operation that may be conceptually simpler and easier to implement experimentally. Of particular interest is to relate our quantum-based model to more classically oriented formulations of CIM, such as those based on coherent state feedback networks with no nonlinearities [Clements], or those based on deterministic nonlinear dynamics (with no quantum noise and only random initial conditions), which have been found to be useful models there to investigate the role of feedback and nonlinearities for CIM combinatorial optimization [45, 47].

[0068] Therefore, MFB-CIMs with non-positive parametric gain can prove useful. Traditionally, MFB-CIMs operate with parametric gain, i.e., pump parameter r>0. However, sampling performance for r≦0 can be investigated, with the r=0 case being of particular experimental interest since it does not require a pump source. Such modifications are easily accommodated within the general Gaussian model, so that the performance can be directly compared against the conventional r>0 case, with gain saturation, quantum noise, etc. held constant.

[0069] Furthermore, ετ nlAn MFB-CIM without nonlinear crystals can be analyzed, which is an MFB-CIM without optical nonlinearities, by setting =0, resulting in a "coherent-state" MFB-CIM since squeezing never occurs. This model has been investigated above in Ref.

[35] in the context of combinatorial optimization (also by a discrete-time formulation), while its performance for Ising sampling can be investigated. Since the resulting system has linear dynamics, the Gaussian formulation applies exactly and is an efficient representation of the quantum state throughout the dynamics.

[0070] Mean-field MFB-CIM with injected measurement noise is a general approach to investigating open dissipative optical systems with weak single-photon nonlinearities, to ignore quantum noise together by adopting the mean-field or classical limit, resulting in a c-number continuous-time differential equation. This limit can also be adopted and motivated for our Gaussian MFB-CIM model, resulting in a new discrete-time model for this mean-field limit, as well as the usual continuous-time mean-field model for MFB-CIM [45, 46]. However, to investigate the sampling performance in this limit, an alternative noise source in the mean-field model is needed. For this, fixed-variance Gaussian distributed noise (which limits to white noise in the continuous-time limit) was injected into the measurement and feedback steps

[48] ; such an extrinsic noise source could correspond, for example, to classical Johnson noise in the detector or random noise intentionally generated by FPGA circuitry (e.g., pseudorandom number generators).

[0071] In FIG. 7(a), the required sampling time as a function of feedback gain is T decay = 4 and η escThe results are shown over a range of pump parameters r, including r≦0 for the above considered case of r=0.2. (For r≧0, these lines are simply vertical slices of the top left panel of FIG. 5.) Surprisingly, the performance is quite comparable over a wide range of pump parameters with negative r corresponding to parametric deamplification, giving slightly better performance at the cost of needing a higher feedback gain to overcome deamplification. In general, for a sufficiently large feedback gain, there is always a robust region over which acceptable sampling performance is obtained, but for r>0, there is also a “sweet spot” at lower feedback gains, where stochastic noise due to antisqueezing can allow efficient sampling with lower feedback gains. In the following, two operation modes, r>0 and r<0, are explored and scaled for larger problem instances.

[0072] A second model is investigated in which the nonlinear crystal is removed from the MFB-CIM, where ετ nl =0 (so, n sat =0), eliminating the need to integrate (16) and (17) for each round trip of the crystal propagation. There are no longer any pump parameters, and T decay and η esc In addition to the feedback gain parameter, we are left with only the feedback gain parameter. Note that in the absence of nonlinear saturation, if the feedback gain exceeds the round-trip attenuation due to losses and light extraction, the system will be unstable, but since our sampling metric (50) only contains the sign of the homodyne result, the metric is unaffected as long as the simulation terminates before numerical overflow. In Figure 7(b), the required sampling time is N decayFor , it is shown that for this linear MFB-CIM, the performance also improves past a certain threshold obtained at small values ​​of the feedback gain compared to the nonlinear MFB-CIM. However, the time required for sampling is at least twice as large as that of the nonlinear MFB-CIM. This observation suggests that nonlinear saturation plays an important role when embedding the Ising problem within the dynamics of the MFB-CIM, and is consistent with the findings of

[47] .

[0073] While the above two cases are straightforward to address within our model, a third approach involves adopting the mean-field limit, which we can motivate as follows. For simplicity, the limit is shown using the continuous-time Gaussian state EOM (27), but by using the exact mapping detailed above, the procedure for the discrete-time version can be derived similarly. To adopt the mean-field limit, the classical mean-field coordinates are

number

number

number

number

number

number

number

number

number

number

[0074]

number

number

number

number

number

number

number

[0075] In FIG. 7(c), the required sampling time is N decay = 4 and η esc = 0.2 as a function of pump parameters and feedback gain for this mean-field model. The left panel shows the mean-field model traditionally used to investigate combinatorial optimization in classical CIM, where sampling is significantly less efficient compared to Gaussian models.

number

number

[0076] Both the coherent-state linear model and the mean-field nonlinear model, each in their own way, explicitly exclude quantum correlations between the internal and out-coupling pulses (i.e.,

number

[0077] The scaling of the sampling performance of the discrete-time MFB-CIM with respect to the problem size can be analyzed and focused on verifying that the different insights and results derived above from investigating small and specific problem instances can be generalized to various as well as large problem instances. The discussion is for a specific Ising problem class with a large number of bases and first excited states, and the performance of the MFB-CIM with multiple instances of this problem class for all given problem sizes can be evaluated. The relationship between the sampling performance and the difficulty associated with a given problem matrix is ​​also explored by investigating how the performance of various alternative models of the MFB-CIM are introduced at a given scale.

[0078] This analysis uses the required sampling time T (used previously) to characterize the sampling performance of the MFB-CIM. samp A more stringent metric can be used because the required sampling time assumes that the machine can be stopped at the optimal time, and while the required sampling time is useful to characterize the possible computational power of the MFB-CIM, it may not represent the actual performance of the MFB-CIM in an experimental setting. all is defined as the number of trajectories required to sample all ground and first excited states multiplied by a fixed number of round trips per trajectory. This definition fits well with the experimental procedure, where each trajectory requires a preselected fixed number of round trips T sim It is simply started for T all gives the time during which an experiment will need to be launched. This choice allows for the use of more sophisticated experiment heuristics to predict when to stop each trajectory, T all It is conservative in the sense that it may result in a faster sampling than given by (and close to the required sampling time). Finally, the time T any is T sumThe number of trajectories required to sample either the ground or first excited configuration, which is similarly defined as the number of trajectories required to sample any of the ground or first excited states multiplied by

[0079] Figure 8 shows how the sampling performance of MFB-CIM scales with problem size N. For any given N, 50 problem instances are considered from the SK1 problem class. A representative instance of this problem class can be found in Figure 8(a), which shows that the off-diagonal elements of the problem matrix are J ij ∈±1. As shown in Fig. 8(b), this problem class has a large total number of degenerate bases and first excitation configurations N conf FIG. 8(c) shows the sampling time T all The distribution of r is shown. Numerical investigations have been carried out on the negative pump (r<0) MFB-CIM model and the model is found to work well for the N=16 instances investigated above. Various parameters associated with the model are specified in Table I.

[0080] [Table 1]

[0081] As shown in the table, some important parameters vary probabilistically from trajectory to trajectory to account for the dynamic threshold variation of the MFB-CIM, leading to some problem instances getting stuck (which seems easy to detect and correct experimentally). Figure 8(c) shows the results of all N confWe also show the wall-clock times required to sample the states. In particular, Fig. 8(c) shows that for a relatively large problem size of N=100, all states can be sampled within 0.6 ms; for future work, it may be interesting to conduct a more thorough benchmark study of how these wall-clock times compare to the wall-clock times achieved by modern algorithms running on conventional digital hardware. Similarly, Fig. 8(d) shows the sampling times T for any basis or first excitation configuration. any The distribution of T any It is noted that T scales exponentially with the problem size N, which is consistent with previous results in [7]. all For scaling of T any Investigating the scaling of,N,provides a better understanding of the overhead incurred when simultaneously sampling,multiple configurations of interest; the difference in exponential base (1.05 vs. 1.08) suggests that the overhead scales exponentially,,but the penalty (with an exponential base of ~1.028) is not particularly high.,Finally, in Figure 8(e), we show that the number of sampling configurations,N, conf With respect to the normalized sampling time

[0082]

number

[0083] Having investigated the details of how the sampling performance of the MFB-CIM scales for a particular operating mode, we consider how different alternative operating modes work with respect to each other. all or any state T anyWe show how the median time to sample r scales with N for the different alternative models

[49] considered above. Due to the large parameter space present in MFB-CIM, the parameters used for each model were heuristically selected by optimizing over a small set of variations (order 10) around the optimal parameters found above for N=16, as listed in Table I. The results confirm that a linear model performs the worst, and thus coherent state measurement noise alone without nonlinearity cannot adequately explain the sampling performance of MFB-CIM. More interestingly, the results also show that setting the pump to be negative or even zero results in sampling performance that scales better for sampling all states as well as any state together. Given the experimental ease of not needing to pump the system, the r=0 results are particularly relevant for future MFB-CIM experiments.

[0084] FIG. 11 shows a computer executing an embodiment of the MFB-CIM system for efficient sampling. Instead of developing quantum computing hardware using superconducting circuits, ion traps, quantum dots, etc. as shown in FIG. 1, the system 1100 uses the computer and theory described above. The models and methods can be programmed into a programmable circuit 1102, which can be an existing digital electronic circuit such as a Field Programmable Gate Array (FPGA) or an Application Specific Integrated Circuit (ASIC). This approach easily launches quantum computing in cyberspace using relatively cheap and reliable silicon integrated circuits rather than requiring costly quantum hardware. In other implementations, the method can be instantiated in a CPU or GPU, and any of these processing cores will execute instructions to perform the method disclosed in FIG. 10.

[0085] In an implementation of the example system shown in FIG. 11, the programmable circuit 1102 to be programmed can be part of a computer system that executes a method on the programmable circuit 1102 to implement the example sampling method discussed above. The computer system can have a display 1104 and a chassis 1106 that interfaces to the programmable circuit 1102 and houses at least one processor 1108 and a memory 1110 from which a number of instruction lines can be executed by the processor 1108 to implement the method stored in the programmable circuit 1102. The chassis can also house the programmable circuit 1102. The example computer system of FIG. 11 can also have typical input / output devices such as a display, a keyboard, and / or a mouse. Alternatively, the computer system can be a server computer, cloud computing resource, etc. that can execute method processes (among other digital storage) on the programmable circuit 1102 to implement the example combinatorial quantum inspired optimizing method discussed below.

[0086] The above description has referred to specific embodiments for purposes of explanation. However, the illustrative discussion above is not intended to be exhaustive or to limit the disclosure to the precise form disclosed. Many modifications and variations are possible in light of the above teachings. The embodiments have been chosen and described in order to best explain the principles of the disclosure and its practical application, so as to enable those skilled in the art to best utilize the disclosure and various embodiments with various modifications as suited to the particular use contemplated.

[0087] The systems and methods disclosed herein may be performed by one or more components, systems, servers, applications, other subcomponents, or may be distributed among such elements. When implemented as a system, such a system may include and / or require components such as software modules, general purpose CPU, RAM, etc., found in a general purpose computer, among others. In embodiments where the innovations reside on a server, such a server may include and / or require components such as a CPU, RAM, etc., found in a general purpose computer.

[0088] Additionally, the systems and methods herein may be accomplished by implementation using disparate or entirely different software, hardware, and / or firmware components beyond those described above. For example, with respect to such other components (e.g., software, processing components, etc.) and / or computer-readable media related to or embodying the present invention, aspects of the innovations herein may be implemented in accordance with numerous general-purpose or special-purpose computing systems or configurations. Various exemplary computing systems, environments, and / or configurations that may be suitable for use with the innovations herein may include, but are not limited to, software or other components embodied in or on personal computers, servers, or server computing devices, such as routing / connectivity components, handheld or laptop devices, multiprocessor systems, microprocessor-based systems, set-top boxes, consumer electronic devices, network PCs, other existing computer platforms, distributed computing environments that include one or more of the above systems or devices.

[0089] In some cases, aspects of the systems and methods may be accomplished or implemented by logic and / or logic instructions, including, for example, program modules executed in conjunction with such components or circuitry. Generally, program modules may include routines, programs, objects, components, data structures that perform particular tasks or execute particular instructions herein. The present invention may similarly be implemented in the context of a distributed software, computer, or circuit configuration, where circuitry is connected by a communication bus, circuitry, or link. In a distributed configuration, control / instructions may originate from both local and remote computer storage media, including memory storage media.

[0090] The software, circuitry, and components herein may also include and / or utilize one or more types of computer-readable media. A computer-readable medium may be any available medium present on, associated with, or accessible by such circuitry and / or computing components. By way of example, and not limitation, computer-readable media may comprise computer storage media and communication media. Computer storage media includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storage of information, such as computer-readable instructions, data structures, program modules, or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory, or other memory technology, CD-ROM, digital versatile disk (DVD), or other optical storage, magnetic tape, magnetic disk storage, or other magnetic storage devices, or any other medium that may be used to store the desired information and that may be accessed by a computing component. Communication media may include computer-readable instructions, data structures, program modules, and / or other components. Additionally, communication media may include wired media, such as a wired network or direct-wired connection, however any such type of media herein does not include transitory media. Combinations of any of the above are also included within the scope of computer-readable media.

[0091] In this description, the terms components, modules, devices, etc. may refer to any type of logical or functional software element, circuit, block, and / or process that may be implemented in various ways. For example, the functions of various circuits and / or blocks may be combined with each other to become any other number of modules. Each module may further be implemented as a software program stored on a tangible memory (e.g., random access memory, read only memory, CD-ROM, hard disk drive, etc.) that is read by a central processing unit to perform the functions of the innovations herein. Alternatively, the modules may include programming instructions transmitted by a communication carrier wave to a general-purpose computer or to processing / graphics hardware. Similarly, the modules may be implemented as hardware logic circuitry that performs the functions encompassed by the innovations herein. Finally, the modules may be implemented using dedicated instructions (SIMD instructions), field programmable logic arrays, or any mixture thereof that provides the desired level of performance and cost.

[0092] As discussed herein, features not inconsistent with the present disclosure may be implemented by computer hardware, software, and / or firmware. For example, the systems and methods disclosed herein may be embodied in various forms, including, for example, databases, digital electronic circuitry, firmware, data processors such as computers also including software, or combinations thereof. Furthermore, while some of the disclosed embodiments describe specific hardware components, systems and methods not inconsistent with the innovations herein may be implemented using any combination of hardware, software, and / or firmware. Furthermore, the features described above, as well as other aspects and principles of the innovations herein, may be implemented in a variety of environments. Such environments and associated applications may include general-purpose computers or computing platforms that may be specifically constructed to perform various routines, processes, and / or operations in accordance with the present invention, or that are selectively activated or reconfigured by code to provide the necessary functionality. The processes discussed herein are not inherently related to any particular computer, network, architecture, environment, or other apparatus, and may be implemented by any suitable combination of hardware, software, and / or firmware. For example, various general-purpose machines may be used with programs written in accordance with the teachings of the invention, or it may prove more convenient to construct a specialized apparatus or system to implement the required methods and techniques.

[0093] Aspects of the methods and systems described herein, such as logic, may similarly be implemented as functions programmed into any of a variety of circuitry, including programmable logic devices (PLDs) such as field programmable gate arrays (FPGAs), programmable array logic (PAL) devices, electronically programmable logic and memory devices, and standard cell-based devices, as well as application specific integrated circuits. Some other possibilities for implementing aspects include memory devices, microcontrollers with memory (such as EEPROM), embedded microprocessors, firmware, software, and the like. Additionally, aspects may be embodied in microprocessors with software-based circuit emulation, discrete logic (sequential and combinatorial), custom devices, fuzzy (neural) logic, quantum devices, and hybrids of any of the above device types. The underlying device technology can be provided in a variety of component types, e.g., metal-oxide semiconductor field-effect transistor (MOSFET) technologies such as complementary metal-oxide semiconductor (CMOS), bipolar technologies such as emitter-coupled logic (ECL), polymer technologies (e.g., silicon-conjugated polymer and metal-conjugated polymer-metal structures), analog-digital mixed, etc.

[0094] It should also be noted that the various logics and / or functions disclosed herein, in terms of their behavior, register transfers, logic components, and / or other characteristics, may be made available using any number of combinations of hardware, firmware, and / or as data and / or instructions embodied in various machine-readable or computer-readable media. Computer-readable media on which such formatted data and / or instructions may be embodied include, but are not limited to, various forms of non-volatile storage media (e.g., optical, magnetic, or semiconductor storage media), but also do not include transitory media. Unless the context clearly requires otherwise, throughout the description, the words "comprise," "comprising," and the like are to be construed in an inclusive sense as opposed to an exclusive or exhaustive sense; i.e., "including, but not limited to." Words using the singular or plural also include the plural or singular, respectively. Additionally, the words "herein," "hereunder," "above," "below," and words of similar import refer to this application as a whole and not to any particular portions of this application. When the word "or" is used in reference to a list of two or more items, the word covers all of the following interpretations of the word: any item of the list, every item of the list, and any combination of items in the list.

[0095] While certain presently preferred embodiments of the present invention have been specifically described herein, it will be apparent to those skilled in the art to which the present invention pertains that changes and modifications of the various embodiments shown and described herein may be made without departing from the spirit and scope of the invention. It is therefore intended that the present invention be limited only to the extent required by applicable legal provisions.

[0096] While the foregoing has been directed to particular embodiments of the present disclosure, it will be recognized by those skilled in the art that modifications of the embodiments may be made without departing from the principles and spirit of the present disclosure, the scope of which is defined by the appended claims. References [1] A. Lucas, Front. Phys. 2, 5 (2014). [2] SG Brush, Rev. Mod. Phys. 39, 883 (1967). [3] F. Barahona, J. Phys. A 15, 3241 (1982). [4] M. Johnson, M. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. Berkley, J. Johansson, P. Bunyk, E. Chapple, C. Enderud, J. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. Thom, E. Tolkacheva, C. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Nature 473, 194 (2011). [5] S. Boixo, TF R nnow, SV Isakov, Z. Wang, D. Wecker, DA Lidar, JM Martinis, and M. Troyer, Nat. Phys. 10, 218 (2014). [6] A. Marandi, Z. Wang, K. Takata, RL Byer, and Y. Yamamoto, Nat. Photon. 8, 937 (2014). [7] McMahon PL, Marandi A, Haribara Y, Hamerly R, Langrock C, Tamate S, Inagaki T, Takesue H, Utsunomiya S, Aihara K, Byer RL, Fejer MM, Mabuchi H, and Yamamoto Y, Science 354, 614 (2016). [8] Inagaki T, Haribara Y, Igarashi K, Sonobe T, Tamate S, Honjo T, Marandi A, McMahon PL, Umeki T, Enbutsu K, Tadanaga O, Takenouchi H, Aihara K,K-i. Science 354 , 603 ( 2016 ). [9] T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, Nat. Photon. Rev. 10, 415 (2016).

[10] I. Mahboob, H. Okamoto, and H. Yamaguchi, Sci. Adv. 2, e1600236 (2016).

[11] T. Wang and J. Roychowdhury, in Unconventional Computation and Natural Computation, edited by I. McQuillan and S. Seki (Springer International Publishing, 2019)pp. 232-2

[12] D. Pierangeli, G. Marcucci, and C. Conti, Phys. Rev.Lett. Rev. 122, 213902 (2019).

[13] Y. Okawachi, M. Yu, J. K. Jang, X. Ji, Y. Zhao, B. Kim, M. Lipson, and A. Gaeta, Nat. Commun. 11, 4119 (2020).

[14] J. Chou, S. Bramhavar, S. Ghosh, and W. Herzog, Sci. Rep. 9, 14786 (2019).

[15] F. Cai, S. Kumar, T. Van Vaerenbergh, X. Sheng, R. Liu, C. Li, Z. Liu, M. Foltin, S. Yu, Q. Xia, J. J. Yang, R. Beausoleil, W. D. Lu, and J. P. Strachan, Nat. Electron. 3, 409 (2020).

[16] Z. Zhu, A. J. Ochoa, and H. G. Katzgraber, Phys. Rev. E 99, 063314 (2019).

[17] G. E. Hinton, Neural Comput. 14, 1771 (2002).

[18] R. Salakhutdinov, A. Mnih, and G. Hinton, in ACM Int. Conf. Proceeding Ser., Vol. 227 (ACM Press, New York, New York, USA, 2007) pp. 791-798.

[19] A. Perdomo-Ortiz, M. Benedetti, J. Realpe-G omez, and R. Biswas, Quantum Sci. Technol. 3, 030502 (2018).

[20] R. S. Bohacek, C. McMartin, and W. C. Guida, Med. Res. Rev. 16, 3 (1996).

[21] V. Lounnas, T. Ritschel, J. Kelder, R. McGuire, R. P. Bywater, and N. Foloppe, Comput. Struct. Biotechnol. J. 5, e201302011 (2013).

[22] K. Ogata, T. Isomura, S. Kawata, H. Yamashita, H. Kubodera, and S. J. Wodak, Molecules 15, 4382 (2010).

[23] H. Sakaguchi, K. Ogata, T. Isomura, S. Utsunomiya, Y. Yamamoto, and K. Aihara, Entropy 18, 365 (2016).

[24] Z. Bian, F. Chudak, R. Israel, B. Lackey, W. G. Macready, and A. Roy, Front. Phys. 2, 56 (2014).

[25] Z. Bian, F. Chudak, R. B. Israel, B. Lackey, W. G. Macready, and A. Roy, Front. ICT 3, 14 (2016).

[26] Y. Matsuda, H. Nishimori, and H. G. Katzgraber, J. Phys. Conf. Ser. 143, 012003 (2009).

[27] S. Mandr a, Z. Zhu, and H. G. Katzgraber, Phys. Rev. Lett. 118, 070502 (2017).

[28] M. S. Konz, G. Mazzola, A. J. Ochoa, H. G. Katzgraber, and M. Troyer, Phys. Rev. A 100, 030303 (2019).

[29] Y. Yamamoto, T. Leleu, S. Ganguli, and H. Mabuchi, Appl. Phys. Lett. 117, 160501 (2020).

[30] S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).

[31] C. M. Caves, Phys. Rev. D 26, 1817 (1982).

[32] H. Wiseman and G. Milburn, Quantum Measurement and Control (Cambridge University Press, 2010).

[33] S. Kako, T. Leleu, Y. Inui, F. Khoyratee, S. Reifenstein, and Y. Yamamoto, Adv. Quantum Technol. 3, 2000045 (2020).

[34] Y. Inui and Y. Yamamoto, Phys. Rev. A 102, 062419 (2020).

[35] W. R. Clements, J. J. Renema, Y. H. Wen, H. M. Chrzanowski, W. S. Kolthammer, and I. A. Walmsley, Phys. Rev. A 96, 043850 (2017).

[36] A. Yamamura, K. Aihara, and Y. Yamamoto, Phys. Rev. A 96, 053834 (2017).

[37] C. Weedbrook, S. Pirandola, R. Garc a-Patr on, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).

[38] G. Adesso, S. Ragy, and A. R. Lee, Open Syst. Inf. Dyn. 21, 1440001 (2014).

[39] J. B. Brask, arXiv:2102.05748 [quant-ph].

[40] I. G. Vladimirov and I. R. Petersen, arXiv:1202.0946 [quant-ph].

[41] C. W. Gardiner and M. J. Collett, Phys. Rev. A 31, 3761 (1985).

[42] H. M. Wiseman and G. J. Milburn, Phys. Rev. A 47, 642 (1993).

[43] S. Lloyd, Science 273, 1073 (1996).

[44] J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics, 2nd ed. (Cambridge University Press, 2017).

[45] Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, Phys. Rev. A 88, 063853 (2013).

[46] T. Leleu, Y. Yamamoto, P. L. McMahon, and K. Aihara, Phys. Rev. Lett. 122, 040607 (2019).

[47] M. C. Strinati, L. Bello, E. G. D. Torre, and A. Pe'er, arXiv:2011.09490 [cond-mat.stat-mech].

[48] D. Pierangeli, G. Marcucci, D. Brunner, and C. Conti, Nanophotonics 9, 4109 (2020).

[49] We do not consider the classical mean- eld model here as realizing those models in experiment requires more energy, making the comparison implicitly unfair.

[50] Y. Inui and Y. Yamamoto, arXiv:2009.10328 [physics.optics]. Appendix A: Orthogonal equations of motion for crystal propagation The total mean field equation of motion is

number

number

number

number

number

number

Claims

1. 1. A system for sampling ground states and low energy Ising configurations using a measurement feedback based coherent Ising machine (MFB-CIM), comprising: a nonlinear optical parametric oscillator configured to be pumped with an internal pulse, the nonlinear optical parametric oscillator using a discrete-time Gaussian state quantum model to generate nonlinear gain saturation of the internal pulse; and a linear measurement feedback loop configured to perform a homodyne measurement of the internal pulse with the nonlinear gain saturation, the linear measurement feedback loop further comprising electronic circuitry configured to perform matrix vector multiplication based on the homodyne measurement.

2. The system of claim 1 , wherein the nonlinear optical parametric oscillator comprises a nonlinear crystal.

3. The system of claim 1 , wherein the electronic circuitry is further configured to calculate a feedback pulse for the nonlinear optical parametric oscillator based on the homodyne measurement.

4. 4. The system of claim 3, wherein an interaction between at least one of the internal pulses and the feedback pulse steers the system toward a lower energy Ising spin configuration.

5. The system of claim 3 , further comprising a coherent injector configured to inject the feedback pulse into the nonlinear optical parametric oscillator.

6. The system of claim 1 , further comprising an outcoupler configured to couple the internal pulse to the linear measurement feedback loop.

7. 2. The system of claim 1, wherein a bistable phase state of the nonlinear optical parametric oscillator due to the pumping of the internal pulse is configured to encode an Ising spin.

8. The system of claim 1 , comprising N nonlinear optical parametric oscillators configured to model an Ising coupling matrix J.

9. 9. The system of claim 8, comprising N linear measurement feedback loops correspondingly coupled to the N nonlinear optical parametric oscillators, each of the N nonlinear optical parametric oscillators configured to generate a corresponding nonlinear gain saturation of a corresponding internal pulse using the discrete-time Gaussian state quantum model.

10. The system of claim 1 , wherein the electronic circuitry comprises a field programmable gate array.

11. 1. A method for sampling ground states and low energy Ising configurations using a measurement feedback based coherent Ising machine (MFB-CIM), comprising: Pumping the internal pulse into a nonlinear optical parametric oscillator; generating nonlinear gain saturation of the internal pulse using a discrete-time Gaussian state quantum model with the nonlinear optical parametric oscillator; performing a homodyne measurement of the internal pulse through the nonlinear gain saturation by a linear measurement feedback loop; and performing a matrix vector multiplication based on the homodyne measurement by electronic circuitry of the linear measurement feedback loop; A method comprising:

12. The method of claim 11 , wherein the nonlinear optical parametric oscillator comprises a nonlinear crystal.

13. The method of claim 11 , further comprising calculating, by the electronic circuitry, a feedback pulse for the nonlinear optical parametric oscillator based on the homodyne measurement.

14. 14. The method of claim 13, wherein an interaction between at least one of the internal pulses and the feedback pulse causes the system to steer towards a lower energy Ising spin configuration.

15. The method of claim 13 , further comprising injecting the feedback pulse into the nonlinear optical parametric oscillator with a coherent injector.

16. The method of claim 11 , further comprising coupling the internal pulse to the linear measurement feedback loop by an outcoupler.

17. 12. The method of claim 11, wherein a bistable phase state of the nonlinear optical parametric oscillator due to the pumping of the internal pulse encodes an Ising spin.

18. The method of claim 11 , wherein the MFB-CIM comprises N nonlinear optical parametric oscillators that model an Ising coupling matrix J.

19. 20. The method of claim 18, wherein the MFB-CIM comprises N linear measurement feedback loops correspondingly coupled to the N nonlinear optical parametric oscillators, each of the N nonlinear optical parametric oscillators generating a corresponding nonlinear gain saturation of a corresponding internal pulse using the discrete-time Gaussian state quantum model.

20. The method of claim 11 , wherein the electronic circuitry comprises a field programmable gate array.

Images