Programmable 200 GOPS Hopfield-inspired photonic Ising machine

TL;DR

The paper addresses the scalability and speed bottlenecks of Ising machines by introducing a Hopfield-inspired photonic Ising machine (CMIM) that operates at room temperature with linear spin representation. It implements the Ising dynamics via a time-multiplexed matrix-vector multiplication in an electro-optic feedback loop, using two cascaded thin-film lithium niobate modulators, a quantum-dot SOA, a controllable optical-noise source for annealing, and a DSP engine to process iterations. The CMIM achieves large-scale problem solving, including up to 256 fully connected spins (65,536 couplings) and 41,209 sparsely connected spins, with >200 GOPS performance and best-in-class solution quality on Max-Cut, protein folding, and number-partitioning benchmarks, outperforming prior photonic Ising machines and rivaling quantum annealers on some tasks while avoiding cryogenic operation and graph-embedding overhead. It demonstrates robust convergence aided by intrinsic noise and DSP-based optimization, and paves the way for scalable, ultrafast optimization, neuromorphic processing, and analog AI using programmable photonics with potential for parallelization and improved energy efficiency.

Abstract

Ising machines offer a compelling approach to addressing NP-hard problems, but physical realizations that are simultaneously scalable, reconfigurable, fast, and stable remain elusive. Quantum annealers, like D-Wave's cryogenic hardware, target combinatorial optimization tasks, but quadratic scaling of qubit requirements with problem size limits their scalability on dense graphs. Here, we introduce a programmable, stable, room-temperature optoelectronic oscillator (OEO)-based Ising machine with linear scaling in spin representation. Inspired by Hopfield networks, our architecture solves fully-connected problems with up to 256 spins (65,536 couplings), and $>$41,000 spins (205,000+ couplings) if sparse. Our system leverages cascaded thin-film lithium niobate modulators, a semiconductor optical amplifier, and a digital signal processing (DSP) engine in a recurrent time-encoded loop, demonstrating potential $>$200 giga-operations per second for spin coupling and nonlinearity. This platform achieves the largest spin configuration in an OEO-based photonic Ising machine, enabled by high intrinsic speed. We experimentally demonstrate best-in-class solution quality for Max-Cut problems of arbitrary graph topologies (2,000 and 20,000 spins) among photonic Ising machines and obtain ground-state solutions for number partitioning and lattice protein folding - benchmarks previously unaddressed by photonic systems. Our system leverages inherent noise from high baud rates to escape local minima and accelerate convergence. Finally, we show that embedding DSP - traditionally used in optical communications - within optical computation enhances convergence and solution quality, opening new frontiers in scalable, ultrafast computing for optimization, neuromorphic processing, and analog AI.

Figures (7)

• Figure 1: Conceptual illustration of the proposed CMIM. Mapping combinatorial optimization problems, such as number partitioning and protein folding, to the Ising model provides an efficient framework for addressing NP-complexity challenges. The minimum energy configurations in the energy landscape ($H$) correspond to solutions to the problem. The hardware comprises two cascaded thin-film lithium niobate (TFLN) modulators lin2024, which encode the input spin vectors and the weight matrix (W). A quantum-dot semiconductor optical amplifier (QD SOA) amplifies the element-wise multiplied signals. A bulk SOA serves as a controlled optical noise source to accelerate convergence when the system struggles to find optimum solutions due to frustrated energy landscapes containing multiple local minima. The optical signal is then converted to an electrical signal using a photodetector (PD). The digital-to-analog converter (DAC) encodes the input signals, while the signal from the PD is sampled by an analog-to-digital converter (ADC), forming an electro-optic (EO) feedback system. Digital signal processing (DSP) stacks process the input and output signals at each iteration. The feedback loop is repeated across multiple iterations until convergence to the ground state solution is achieved. The mathematical representation (top) shows how the weight matrix is determined from the fixed coupling matrix $\mathbf{J}$ and the local bias $\mathbf{h}$ corresponding to the given problem. The coupling matrix is flattened for element-wise multiplication, and the resulting values are summed to perform matrix-vector multiplication. The DSP processes the received data and routes it back through the feedback path. Over successive iterations, the system converges to the ground state solution ($\sigma$).
• Figure 1: Detailed experimental schematic of the CMIM.. The electrical signals were generated using a two-channel Keysight M8199B arbitrary waveform generator (AWG) as the DAC with a sampling rate of 256 GSa/s. Channel 1 encoded the Ising spins, while Channel 2 encoded the fixed weight matrix. The TFLN modulators have a measured 3-dB bandwidth of 110 GHz and a DC half-wave voltage of $V_\pi=1.5$ V. The optimized transmitter (Tx) and receiver (Rx) digital signal processing (DSP) stack have been employed during high-speed operation. The modulated optical signal was detected using a PD with a 3dB bandwidth of 100 GHz, and sampled by a 256 GSa/s Keysight UXR real-time oscilloscope (RTO) as the ADC.
• Figure 2: Experimental characterization of bifurcation and matrix-vector multiplication. (a) Heatmap of the bifurcation dynamics as a function of feedback strength while operating at 64 GBaud with $N=$ 262,144 uncoupled spins. The colour indicates the number of spins in a given state after 50 iterations. (b) Heatmap showing the evolution of uncoupled spins over successive iterations at a feedback strength of $\alpha=3.5$ (left), alongside a histogram of spin values at the final iteration (right). (c) Experimental results of matrix multiplication performed at 64 GBaud using 500 randomly sampled 128$\times$128 matrices. (d) Bit precision of matrix multiplication as a function of matrix size, ranging from $N=32$ to $N=200$. (e) As the symbol rate increases from 4 GBaud to 148 GBaud, the bit precision of matrix multiplication decreases linearly from 5.0342 to 2.7885.
• Figure 2: Optical noise injection with a bulk SOA. (a) ASE output power directly out of the bulk SOA and noise variance after measuring with noise with an RTO according to the schematic shown in Fig. \ref{['fig:expSetup']}. (b) Fit Gaussian probability density functions of the measured noise at the RTO for different SOA bias currents. (c) Measured OSNR at varying bulk SOA bias currents; waveform is data modulated for a 400-node lattice problem. The inset shows the wavelength spectrum for two different SOA biases. (d) Iterations required to reach ground state for a 400-node 2D lattice problem, depending on the bulk SOA bias current running at 64 Gbaud. (e) Evolution of the energy of a 400-node 2D lattice problem with constant noise (0 mA and 60 mA bulk SOA bias) and exponential annealing (initial SOA bias is 100 mA).
• Figure 3: Benchmarking results. Experimental results for the square lattice, G22 graph, and G81 graph problems, all evaluated at the maximum operating baud rate of 106 GBaud. (a) Time evolution of the Ising energy during ground state search for a 2D square lattice with 10,201 spins ($101\times 101$ lattice). The inset shows snapshots of the spin domain evolution across iterations, illustrating convergence towards the ground state solution. (b) Scaling analysis of solution quality--- measured by proximity to the ground state---across problem sizes ranging from 100 to 41,209 spins ($203\times 203$ lattice length). (c) Time evolution of the cut value for the G22 graph. The system achieves 99.48% of the best-known cut, reaching a cut value of 13,289. The inset shows the corresponding coupling matrix with 1% edge density and binary weights of 0 and 1. (d) Time evolution of the cut value for the G81 graph. The system reaches a solution quality exceeding 96.34%, with a cut value of 13,516. The inset shows the G81 planar graph with weight values of $+1$ and $-1$. The evolution of cut values for both the G22 and G81 graphs is based on a fixed iteration count of 800.