Entropy Computing: A Paradigm for Optimization in an Open Quantum System

Abstract

Modern quantum technologies using matter are designed as closed quantum systems to isolate them from interactions with the environment. This design paradigm greatly constrains the scalability and limits practical implementation of such systems. Here, we introduce a novel computing paradigm, entropy computing, that works by conditioning a quantum reservoir thereby enabling the stabilization of a ground state. In this work, we experimentally demonstrate the feasibility of entropy computing by building a hybrid photonic-electronic computer that uses measurement-based feedback to solve non-convex optimization problems. The system functions by using temporal photonic modes to create qudits in order to encode probability amplitudes in the time-frequency degree of freedom of a photon. This scheme, when coupled with electronic interconnects, allows us to encode an arbitrary Hamiltonian into the system and solve non-convex continuous variables and combinatorial optimization problems. We show that the proposed entropy computing paradigm can act as a scalable and versatile platform for tackling a large range of NP-hard optimization problems.

1 Introduction

In the domain of computational optimization, a significant number of problems are NP-hard, which is commonly considered the hardest class of optimization problems. It has been shown that there exists a polynomial-time mapping between any problem in NP to an NP-hard problem [1, 2]. It is strongly suspected, but not formally proven, that no efficient (polynomial time) algorithm exists for NP-hard problems, even for quantum computers. However, even if such problems are likely impractical to solve at very large scales, finding high quality approximations (or optimal solutions) in practical timescales is very important, since numerous real world applications can be mapped onto NP-hard problems. The inherently complex nature of such optimization problems necessitates the development of novel computational frameworks, both in terms of hardware and algorithms. Traditional computational methods implemented on classical Complementary Metal Oxide Semiconductor (CMOS) devices often struggle with the scale and complexity of these problems, leading to extended solution times or, in some cases, the inability to find an optimal solution within a reasonable time frame.

Various special purpose processors (both classical and quantum), such as analog solvers, have recently been proposed. A detailed review of some of these machines can be found in [3]. Many of the recently proposed classical special purpose processors are a class of CMOS devices which effectively implement simulated annealing and related algorithms, a review of these devices can be found here [4, 5]. Similarly, many quantum devises used for non-convex optimization are quantum annealers, which were first proposed by Kadowaki and Nishimori [6]. Quantum annealers differ from the more conventional gate-based quantum computing architectures, since they encode and solve problems via continuous time evolution. Recent theoretical investigations showed annealing-based architectures might give equivalent advantages to those seen in a gate model setting. Adiabatic quantum computing, an idealized setting of quantum annealing where very long anneal times are employed, was shown to give the same quadratic speedup [7] that Grover search yields in a gate-model device [8, 9], and is the best that is achievable by any quantum algorithm [10]. Furthermore, a more general class of continuous-time algorithms can obtain the same speedup, including continuous time quantum walks [11] and interpolations between adiabatic and quantum walk protocols [12]. Even more recently, the same advantage has also been shown to be possible through dissipative dynamics, similar to those that our device uses [13]; similar effects can be used to solve optimization problems [14]. The physical realization of quantum annealers on superconducting platforms remains a challenge due to limitations on connectivity and scalability.

A popular quantum analog approach from the past decade, which bears some architectural similarities to the devices we study here, is the setting of coherent Ising machines (CIM) [15, 16, 17] due to its energy efficiency, scalability, and global connectivity. Superior time-to-convergence vs. problem size scaling compared to annealers was observed [18]. However, maintaining stable operation on CIM to avoid external perturbations, including amplitude heterogeneity or phase-stability demand over long fiber, is the main drawback that prevents this technology from being widely adopted [19, 20].

It is important to note that the aforementioned efforts solely focus on solving Ising-type problems. However, the landscape of computationally-hard optimization problems is not limited to Ising problems and contains, for example, binary, non-convex continuous-variable, integer, and mixed-integer problems. Many discrete NP-problems do not naturally or directly map to the Ising model’s framework of binary spin states with quadratic interactions, and developing an effective mapping that accurately represents the original problem within the Ising framework can be highly non-trivial and is often achieved only with a high computational overhead, which renders the approach impractical. Furthermore, incorporating constraints, which are often required in NP-hard problems, into the Ising model can increase the complexity of the problem formulation. This is typically accomplished by introducing additional spins (qubits in the context of quantum computing) and carefully designing the interaction terms to ensure that the constraints are properly enforced, which can significantly increase the size and complexity of the resulting Ising problem. An additional complication, particularly for analog computers, is that the ratio of the largest to the smallest non-zero elements (the dynamic range) can be large in some problems, particularly when constraints are considered, or where mapping procedures to an incompatible hardware graph, for example through minor embedding, must be performed [21]. The small dynamic range of any real hardware often limits the types of problems that can be efficiently solved on that hardware.

In light of these challenges, we introduce a new computing paradigm, entropy computing, that operates by conditioning a quantum reservoir to promote the stabilization of the ground state in a quantum optical setting. In this approach, a target Hamiltonian is mapped onto an effective dissipative operator that hosts the propagation of photon qudit states to solve an optimization problem in the photon number Hilbert space. This approach induces loss or decoherence into the system to suppress the evolution of unwanted states while promoting the evolution of qudits that represent lower energy states of the corresponding Hamiltonian. We experimentally demonstrate entropy computing through a hybrid optoelectronic measurement-feedback system that utilizes photon qudits encoded as probability amplitudes in time-frequency degree of freedom in conjunction with electronic interconnects for embedding an arbitrary Hamiltonian. In this manner, we demonstrate an entropy computing machine with a versatile polynomial loss function containing first- to fifth-order terms with fully programmable weight tensors, capable of solving optimization problems with up to 949 independent variables over a fixed summation constraint. We employ this platform for solving non-convex continuous variables and integer combinatorial optimization problems.

2 Results

2.1 All-optical Entropy Computing System

The working principle of the system is schematically depicted in Fig. 1(a). It consists of an optical-feedback loop that includes an optical amplifier, a photon-mode mixer/encoder, and a loss medium to implement both linear and nonlinear loss mechanisms. Quantum states are encoded as photon numbers in a train of time-bin optical pulses in the loop, which could be either in fiber or free space. During each loop, optical signals are amplified and sent into the mixer that performs linear transformation of the time-bin modes according to the desirable Hamiltonian $\hat{H}$. An example of the mixer uses a series of beamsplitters, optical delay lines, and optical switches, all controlled opto-electronically to construct a multiport, delay-switched Mach-Zehnder interferometer. The mixed signal is then sent to the loss medium, so that the optical system realizes the imaginary time evolution of an open quantum system with linear and nonlinear losses. This loss medium can be implemented via Quantum Zeno blockade where the loss must be stronger than the cavity coupling[22, 23, 24], or sum frequency generation, by which the system would gradually relax to the effective ground state. Let $\{\ket{\Psi_{n}}\}$ be the set of eigenstates of the Hamiltonian as encoded into the system by the mixer and losses, with their corresponding eigen-energies being $E_{n}$. Consider the system to initially be in a mixture state, described by a density matrix of

$$\rho(t=0)=\sum_{n}P_{n}\ket{\Psi_{n}}\bra{\Psi_{n}},$$ (1)

where $\{P_{n}\}$ are the normalized probabilities. Because the system evolves at an imaginary time, the state at a later time $t$ is

$$\rho(t) =\frac{1}{\it{N(t)}}e^{-i\hat{H}(-it)}\rho(0)e^{iH(it)}$$ (2)

$$=\frac{1}{\it{N(t)}}\sum_{n}P_{n}e^{-2E_{n}t}\ket{\Psi_{n}}\bra{\Psi_{n}},$$

where $N(t)$ is the normalization factor. Note that the above equation is an effective description of the open quantum system in the Markovian limit.

Figure 1: An example realization of entropy quantum computing using time-energy modes and its schematic in a fiber-optical loop. (b) The working principle. It is based on carving a desired Hamiltonian ($\hat{H}$) into a dissipation operator ($-\mathrm{i}\hat{H}$) which tends to favorably support the evolution of a quantum state of light that resembles the lower-energy states, particularly the ground state, of the target Hamiltonian, while simultaneously undermining other states through decay/decoherence. (a) The time-domain evolution can be discretized in a fiber-optical loop setting to iteratively evolve the state and relax the system. The quantum states are encoded into a train of time-bin states of light in the photon-number Hilbert space. During each loop, optical signals are coupled with vacuum fluctuation and amplified by an optical amplifier with fixed average power. The output is sent into the mixer and encoder that performs linear transformation of the time-bin modes according to the desirable Hamiltonian. The mixed signals are sent to the loss medium—which causes differential loss to each time mode—before subject to the constant-power amplification. (c) Evolution of the state vector during optimization. This figure illustrates the state vector’s behavior as the system iteratively evolves over time. States with higher loss gradually diminish, while states with lower loss persist and are amplified.

2.2 Hybrid Entropy Computing System

In the first realization of this hardware, which is discussed here, we use time correlated single photon counting (TCSPC) and electro-optical feedback to emulate entropy computing, although for discounted computing power. The system configuration is illustrated in Fig. 2, where the loss mechanism is implemented through an electro-optical modulator (EOM) and the mixer is realized by passing the optical signals through a nonlinear optical circuit, detecting them in a single photon detector, and post-processing the TCSPC results. The optical signal is generated by a continuous-wave laser followed by a variable optical attenuator. It passes through the EOM, the nonlinear circuit, and detected by a single photon detector. The TCSPC results are fed to a field-programmable gate array (FPGA) board, where the radio-wave input to the EOM is calculated and generated in a digital-to-analog converter (DAC) for each time bin. Here, the nonlinear circuit, photon detector, and the FPGA function together as a mixer/encoder.

In this system, the temporal bins form the state bases and TCSPC measures photon counts in each bin for stochastic computing with integer variables [25, 26]. Rather than the amplification of spontaneous emission as in Fig. 1, the quantum fluctuation comes from the shot noise of single photon counting. Through the feedback, the photon counts of each time bin induce subsequent conditional photon losses in the same or other time bins. Hence, the quantum states evolve step-wise through a measurement-and-feedback setting. By enforcing a normalization condition for the total photons in all time bins, the photon numbers in each time bin will be incentivized to self-align into a configuration with the minimum loss, as the high loss configurations will be “punished” through the loop iteration. As we can conveniently control the loss for individual time bin and the overall loss, one can apply constraints to photon numbers in each time bin or the sum of all. This gives the flexibility of optimizing practical problems that always come with various constraints. Throughout this process, the injection of quantum fluctuations enables a search in a near-infinite search space, which assists in bypassing the trapping of the quantum state into local minima. This emulation system, however, is not as efficient as the full system in Fig. 1 in searching the whole Hilbert space. It requires many incremental loops to find the global minimum solutions to optimization problems with many variables.

Figure 2: An emulation system for entropy computing using time-energy modes and a measurement-feedback scheme. (a) The Hamiltonian problem is encoded into the amplitude of an electrical signal via digital-to-analog converter. The signal is used to drive an EOM device which tailors a weak coherent state into single-photon signals in a shaped wave function, to realize high dimensional time-bin encoding. The optical signal is then combined with and modulated by a coherent pump at a different wavelength as they pass through a quantum non-linear optical circuit. The modulation includes but is not limited to sum-frequency generation, quantum Zeno blockade [22, 23, 24], quantum optical arbitrary waveform generation [27, 28]. At the output, the resultant single photons are detected by a single photon detector (SPD) and recorded by a time-to-digital converter (TDC). A clock signal is used as reference for the TDC where the period matches with the feedback loop time. Field-programmable gate array (FPGA) accumulates photon counts of each time-bin and computes the contribution of those counts to the losses of each time bin, thereby emulating the interaction terms in the Hamiltonian. Hence, the loss rate for each mode is the sum of a constant term, corresponding to the linear “chemical potential energy” of that mode, and a photon-number dependent term, corresponding to the nonlinear interaction energy. The quantum non-linear optical circuit and time-correlated single photon counting by the TDC, SPD and feedback through the EOM together act as the medium with linear and nonlinear losses. Single photon detection combined with FPGA normalization feedback promotes least loss states in a hybrid framework. A secondary digital-to-analog converter embedded in the FPGA is used to monitor and control the mean photon number to guarantee high probability of single photon level as the photon wave-function changes every feedback loop. (b) The quantum states are encoded into a train of time-bin states of light in the photon-number Hilbert space. The linear loss in the Hamiltonian is mapped into probability amplitudes of a wave function. Each variable of the objective function is assigned into each time bin, creating “qudit” equivalent that is widely used in high-dimensional temporal encoding in quantum random number generation, quantum key distribution, and single-photon sensing [29, 30, 31]. The probability amplitude of each time bin represents its linear loss coefficient $C_{i}$.(c) Single photon collapsed state is collected one by one and accumulate to create the next feedback to tailor new wave functions in each loop. When the wave function evolutions start freezing, number of photons collected in each bin are translated directly as a state vector multiplied with a normalization factor.

2.3 Dirac-3 Quantum Optimization Machine

The emulation system can be used for a variety of optimization tasks, including those of binary variables, mixed-integer variables, quasi-continuous, and any combination of them. In this work, we focus on the minimization of the following cost function $E$ over variables $v_{i}$:

$$E=\sum_{i}C_{i}v_{i}+\sum_{i,j}J_{ij}v_{i}v_{j}+\sum_{i,j,k}T_{ijk}v_{i}v_{j}v_{k}+$$

$$\sum_{i,j,k,l}Q_{ijkl}v_{i}v_{j}v_{k}v_{l}+\sum_{i,j,k,l,m}P_{ijklm}v_{i}v_{j}v_{k}v_{l}v_{m}~{}.$$ (3)

Here, $v_{i}$ ($i=1,2,3,\cdots,N$) are real numbers over a discrete space, $C_{i}$ are the linear terms’ real-valued coefficients while $J_{ij}$, $T_{ijk}$, $Q_{ijkl}$, $P_{ijklm}$ represent two-body to fifth-body interaction coefficients that are real numbers subject to the tensors $J$, $T$, $Q$, and $P$ being symmetric under all permutations of the indices. Furthermore, it is important to note that the probabilistic nature of the optimization variables in the proposed entropy computing method demands that all variables are non-negative, i.e., $v_{i}\geq 0$. To allow negative variables, a linear transformation of the variables needs to be applied prior to the problem encoding.

In observation of the system openness due to effective dissipation and gain applied during each feedback loop, a constraint on the sum of the variables $v_{i}$ can be applied such that:

$$\sum_{i}v_{i}=R.$$ (4)

Thus, the proposed scheme finds minimal solutions of the loss function of Eq. (3) subject to the sum constraint of Eq. (4).

Given a desired optimization problem formulated in the form of relation (3), the optimization process is described below. First, a sum constraint $R$ is chosen. Clearly, without a knowledge of the final solution, the true sum of variables $\sum_{i}v_{i}$ is unknown. But, one can readily bypass this problem by introducing slack variables, as discussed in an example in the next section. Once the sum is predetermined, an initial photon qudit state is generated from the noise and launched in the system. After propagating in the closed loop, the photonic state is measured and evaluated to prepare the state for the next iteration.

In stark contrast with the Ising Hamiltonian, which is the basis model for the majority of quantum annealers, the above objective function involves polynomial terms (up to fifth order in our present experimental demonstration) over discretized variables. In this regard, the proposed hybrid quantum optimization machine offers two immediate advantages over an Ising solver; (i) it can naturally represent higher than binary optimization problems, and (ii) it involves k-body interaction terms ($k=1,2,3,4,5$). Accordingly, it offers great potential in efficiently solving continuous and integer variables as well as problems that naturally involve higher-order interaction terms such as the satisfiability boolean, without requiring additional complex encoding or incorporating auxiliary variables that add to the size of the problem in case of an Ising solver. Furthermore, the proposed machine naturally allows for dense long-range interactions in all orders of the interactions, which alleviates the requirement for complex embedding algorithms. Here, we report results from our first commercially available machine, which we call Dirac-3. Dirac-3 is an optimization solver which implements the entropy computing paradigm discussed above.

2.4 Experimental Results

In the following, we first use a simple 2-variable problem to illustrate the optimization process in the proposed entropy computing paradigm. Next, we present results from two benchmark hard optimization problems from known libraries. In general, there is no simple method for benchmarking the proposed computing method and system. Here, we use two different types of problems, a non-convex quadratic problem, and a combinatorial graph cut problem.

We first consider a simple two-variable optimization problem, defined through the objective function $g(x,y)=(x^{4}+y^{4})/4-5(x^{3}+y^{3})/3+3(x^{2}+y^{2})$. As depicted in Fig. 3(a), this function contains three local minima at $(\textup{x},\textup{y})=\{(0,3),(3,0),(3,3)\}$ as well as a global minimum at $(\textup{x},\textup{y})=(0,0)$. To solve this problem using Dirac-3, we consider 3 variables $v_{1}$, $v_{2}$, and $v_{3}$, with $v_{1}$ and $v_{2}$ representing $x$ and $y$ respectively, while $v_{3}$ representing a slack variable which is not coupled to any other variable but used to effectively relax the sum constraint. In this manner, one can consider a relatively large sum constraint of $R=100$ to incorporate a larger range of potential solutions. Figure 3(b) shows the evolution of the energy level for 20 different runs as the system evolves toward an equilibrium point. The evolution of the variables versus the number of iterations are shown in Figs. 3(c) for three exemplary cases. These figures indicate rapid convergence of the solution after only a few iterations. Even more interestingly, the trajectories of the optimization variables in the two-dimensional $(x,y)$ plane (Figs. 3(d-k)) show that the system takes steps to avoid the local minima, even though in all cases the initial point was located within the attractor basins of the local minima.

Figure 3: Solving a two-variable non-convex quadratic optimization problem. A two-variable non-convex polynomial optimization problem is considered. (a) A visualization of the energy landscape that involves three local minima and a global minimum at $(\textup{x},\textup{y})=(0,0)$. (b) The iterative evolution of the cost function of the proposed hybrid entropy computing solver over 20 runs. (c) Three exemplary evolutions of the optimization variables involved, including the slack variable ($x_{3}$), over iterations of the entropy computing solver. (d-k) Eight exemplary trajectories of the optimization variables in the two-dimensional $(x,y)$ plane as the solver evolves toward equilibrium while starting from different initial conditions. In these figures, the solid lines show the trajectory of the entropy-computing solver, while the dashed lines depict the trajectory of a gradient-descent solver.

We next consider a non-convex quadratic optimization problem (QPLIB$\_$0018) with 50 continuous variables over a fully connected weighted graph, selected from QPLIB, a library of quadratic programming instances [32]. The cost function of this problem is of the form $f(\mathbf{x})=C^{T}\mathbf{x}+\mathbf{x}^{T}J\textbf{x}$, where, $\textbf{x}=(x_{1},x_{2},\cdots,x_{50})^{T}$ and $x_{i}\geq 0$. Figures 4(a,b) show the equilibrium energy distribution over 50 runs of Dirac-3 compared with the results obtained from a conventional gradient descent algorithm. Here, linear terms are added to the problem, which produces an offset of 1 from the original formulation’s objective value for all solutions. Dirac-3 successfully finds the ground state in $84\%$ of instances. The evolution of the energy level is plotted as a function of the iterations (red) and compared with those obtained from gradient descent (blue) over 50 different runs, shown in Fig. 4(c). The inset depicts the evolution of the solution with the number of iterations indicating rapid convergence of the solution after a few iterations.

Figure 4: Solving a non-convex continuous optimization problem. A non-convex quadratic optimization problem with 50 variables is considered (QPLIB$\_$0018). (a,b) Energy distribution over 50 runs of Dirac-3 (red) and gradient descent algorithm (blue). (c) Energy evolution versus the number of iterations on Dirac-3 (red) and gradient descent (blue). The inset shows the evolution of the solution versus iterations. (d), (e) Investigate the relationship between the mean photon number ($\mu$) and two key performance metrics: returned energy and time-to-convergence. The blue gradient, transitioning from light to dark, represents a decreasing mean photon number ($\mu$). This signifies that Dirac-3 is operating progressively closer to the single-photon regime.

It is important to evaluate the performance of the entropy computing machine while operating with classical versus quantum states of light For this purpose, we investigate the optimization performance for various mean photon numbers $\mu_{p}$ associated with both regimes. To do so, in each feedback loop, we accumulate single-photon counts. Here, we experiment with different mean photon numbers per run, which consists of multiple feedback loops. A lower mean photon number signifies low probability of multi-photon in each coherent light pulse. This concept has been used widely in coherent state quantum key distribution where, if $\mu_{p}$ is too large, information is encoded redundantly in multiple photons, which damages the security of secret key exchange [33, 34]. In this case, a high probability of operating in the single-photon regime is required in every feedback loop to constantly boost the quantum stochasticity of the system, promoting chances of jumping out of local minima. We varied $\mu_{p}$ from $3$dB below to $3$dB above its experimental optimal level. Figure 4(d) shows the equilibrium energy level versus the mean photon number. This figure indicates the higher probability of finding the ground state solution for an optimal range associated with low mean photon numbers. By increasing $\mu_{p}$ the solution quality deteriorates as the system is more likely to be trapped in local minima. On the other hand, for lower $\mu_{p}$ the probability of collecting a single photon in every feedback loop increases, which leads to higher probabilities of escaping local minima and converging to the global minimum. Evidently, the evolution trajectory of Dirac-3 in searching for optimal solution in a non-convex energy landscape shows “quantum-tunneling-like” behavior in Fig.  4. Notice that further reducing the single photon rate the classical (thermal) noise from the single-photon detector’s dark count deteriorates the results, too. This can be quantified into a parameter, the quantum-to-classical ratio (QCR), defined as the photon detection rate divided by the dark count rate, which decreases for lower photon numbers. Furthermore, it is important to note that there is a trade-off between increasing the quality of the solution by operating in the quantum regime and reducing the optimization time. As shown in Fig. 4(e) operating at lower photon numbers requires a longer time to accumulate enough photons to overcome high Poisson noise, leading to an increase in time for the system to stabilize in an equilibrium state.

Another critical parameter is the number of accumulated photons (or the number of feedback loops) before updating a new wave function (new iteration). This is directly related to the amount of dissipation in this open quantum system set up. For practical purposes, Dirac-3 is designed to be configurable with different mean photon numbers and different accumulated photons per iteration, which assembled into four relaxation schedules simply named 1,2,3, and 4. As a result, the probability of finding an optimal solution is much higher in schedule 4 than 1, especially on a competitive energy landscape.

Coupling terms that represent nonlinear loss can be realized on MZI mesh, but it is not scalable for large problem size until photonic integrated circuits are available. In this hybrid quantum optimization machine implementation, FPGA can be conveniently engineered to accommodate these higher order of interaction terms, making this architecture more practical with current technology [35, 36]. Thus, on this hybrid measurement-feedback based approach, time-to-convergence depends partially on the FPGA and its architect for multiplier.

Next, we consider a combinatorial optimization problem, the maximum cut (max-cut) problem, that is known to be NP-hard. Considering for simplicity an unweighted graph, the max-cut is a partition of its vertices into two complementary sets, such that the number of edges between the two sets is maximal. This problem can be generalized to maximum k-cut (max-k-cut), where, the vertices are partitioned into $k$ disjoint subsets, such that the number of edges between the $k$ subsets is maximized. Although the proposed entropy computing paradigm deals with continuous-valued variables in general, it can be utilized to solve such discrete problems by suitable encoding of a discrete (e.g., spin) degree of freedom in the continuum, which can be done in different ways. One way is to map a binary variable $s_{i}\in\{0,1\}$ onto a vector $\vec{\mathbf{s}}_{i}=x_{i}\bigl{(}\begin{smallmatrix}1\\ 0\end{smallmatrix}\bigr{)}+y_{i}\bigl{(}\begin{smallmatrix}0\\ 1\end{smallmatrix}\bigr{)}$. In this manner, one can embed a classical bit onto the continuum and identify the two states by taking the maximum of $x_{i}$ and $y_{i}$ variables. Furthermore, the values of $x_{i}$ and $y_{i}$ can be regularized by considering a regularization term that can be chosen as $x_{i}+y_{i}=1$ in its simplest way. This approach can be readily extended to consider a ternary digit (trit), a quaternary digit (quit), and so on, which allow for tackling a general $k$-state standard Potts problem with the proposed entropy computing machine.

Now, considering a graph with $N$ nodes, described with adjacency matrix $A$, i.e., $A_{ij}=1$ for adjacent nodes and $A_{ij}=0$ otherwise, the max-k-cut problem can be formulated as minimizing the following objective function: $E=\sum_{i,j}A_{ij}\vec{\mathbf{s}}_{i}\cdot\vec{\mathbf{s}}_{j}+\lambda\sum_{i}{\left\|\vec{\mathbf{s}}_{i}-\mathbf{1}\right\|}_{1}^{2}$, where the second term is an L1-norm regularization term with a parameter $\lambda$. It is straightforward to cast this objective function in the form of relation 3. This results in a quadratic function involving $qN$ variable, where $q$ is the dimension of a single artificial spin state and $N$ is the size of the graph.

To test this formulation, we consider a randomly generated graph of $N=30$ nodes such that each two nodes are connected with a probability of $p=0.5$. Figure 5(a-c) visualizes exemplary partitioning of this graph into 2, 3 and 4 groups through the proposed entropy computing system. This particular graph has $233$ total links, and we found the maximum cut size to be $146$. Histograms of the cut sizes obtained for 100 different runs of the max-cut, max-3-cut, and max-4-cut problems using the entropy computing solver at schedule 1 are shown in Figs. 5(d-f), and using relaxation schedule 4 in Figs. 5(g-i). The results are compared with those obtained from Semi-Definite Programming (SDP) shown in Figs. 5(j-l). As expected, we observe that schedule 4 provides better results. In addition, as this figure clearly indicates, in all cases the entropy computing machine provides excellent results that outperform SDP results. As a reference, it is worth mentioning that SDP relaxation for max-cut has a performance guarantee of the ratio of $0.878$ [37], which translates to a cut size of $128$ for the particular graph considered here. For the results presented in Fig. 5, the regularization parameter was chosen to be $\lambda=5$. In general, compared to using a penalty-based approach, we find that the results are less dependent on the regularization parameter and the sum constraint.

Figure 5: The results for solving the max-k-cut problem. The optimization results for solving combinatorial problems of max-cut, max-3-cut, and max-4-cut on a 30-node unweighted graph that is generated randomly with $p=0.5$ probability of connectivity between each two nodes. (a-c) The visualizations of exemplary graph cuts with different numbers of partitions. (d-f) The distribution of the cut sizes over 100 runs of the max-cut, max-3-cut, and max-4-cut problems on the EQC system using relaxation schedule 1. (g-i) Similar distributions when using the schedule 4. (j-l) The distribution of the cut sizes over 100 runs of the same problems using semi-definite programming (SDP).

3 Discussion and Conclusion

We introduced the concept of entropy quantum computing as a novel computing paradigm for solving hard optimization problems. We experimentally demonstrated entropy computing through a hybrid photonic-electronic system that builds on time-multiplexed photon qudit time bins propagating in a closed feedback loop involving electronic interconnects for implementing a reconfigurable effective optimization cost function. We demonstrated the successful operation of the proposed system for solving non-convex and combinatorial optimization problems. Results show that Dirac-3 found the ground state more often than classical gradient descent on a non-convex optimization problem with constraint. Dirac-3 also performs superior to semi-definite programming in solution quality on Ising and standard Potts problems. The presence of the sum constraint in the current Dirac-3 quantum machine becomes an advantage in problems that intrinsically involve this natural restriction. Such problems emerge in portfolio optimization, resource allocation problems, diet problems, knapsack, network flow problems, and election and voting systems [38, 39, 40, 41]. One of the key strengths of the proposed entropy computing machine is its flexibility in encoding, allowing it to handle continuous and integer variables. This capability sets it apart from many classical and quantum solvers that are primarily designed for binary Ising/QUBO problems, and it opens up new possibilities for efficient solutions in diverse applications, including grid optimization and machine learning tasks like clustering and decomposition [42, 43]. Traditionally, implementing higher-order interactions in Ising, Potts, or XY problems on analog hardware solvers requires a quadratization step for polynomial reduction [44, 45, 46]. Dirac-3 simplifies this process by directly mapping high-order optimization problems, promising increase in precision and solution quality while consuming fewer resources by eliminating the need for auxiliary variables [47, 48]. Further studies and benchmarking are necessary to fully explore the capabilities of Dirac-3.

The proposed entropy computing system harnesses noise to steer clear of local minima in complex energy landscapes. Our findings suggest that operating near the single-photon regime (what we call the “quantum” regime) improves the optimization results as the mean photon number decreases. Effectively, time-to-convergence increases as fewer photons require more samples to be collected, since the detection rate will be lower. Faster modulation and higher single photon detection efficiency are key avenues to improve overall computation speed. The current hybrid architecture of the entropy computing system exhibits low energy consumption, comparable to a laptop. Further benchmarking needs to be done to appreciate energy advantages of this architecture. This energy footprint is expected to decrease significantly when the system is implemented entirely on an integrated photonic platform. This shift towards an all-optical approach holds promise for a practical and sustainable unconventional computing paradigm, contributing to solving energy rising issues by high performance computing [49, 50].

Acknowledgements

Authors would like to thank other scientists at Quantum Computing Inc including Prajnesh Kumar, Jeevanandha Ramanathan, and Malvika Garikapati for their valuable engineering suggestions. NC was fully supported by QCi in performing this work, in particular no UKRI support was received for writing this paper.

Declarations

Code Availability. The codes that support the findings of this study are available from the corresponding author upon reasonable request.

Data Availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author Contribution. All authors significantly contributed in preparing the results presented in this manuscript.