PaperlayerView exact revision on arXiv

Source-linked AI summary

Heuristic Recurrent Algorithms for Photonic Ising Machines

Charles Roques-Carmes, Yichen Shen, Cristian Zanoci, Mihika Prabhu, Fadi Atieh, Li Jing, Tena Dubcek, Chenkai Mao, Miles R. Johnson, Vladimir Ceperic, John D. Joannopoulos, Dirk Englund, Marin Soljacic

arXiv:1811.02705v3physics.app-phcs.ETphysics.optics

TL;DR

Heuristic methods struggle with large combinatorial Ising problems, motivating algorithms suited to parallel hardware. PRIS uses recurrent, noisy updates with eigenvalue dropout; the authors report high-probability ground-state finding and use it to probe critical behavior.

  • Problem

    Heuristic methods struggle to solve large problems, while finding the minimizing spin distribution for an arbitrary Ising coupling matrix is NP-hard.

    Source: INTRODUCTION · p. 2

  • Method

    PRIS is a photonic-oriented heuristic that recurrently transforms Ising spin states using matrix operations, dynamic noise, and eigenvalue dropout.

    Source: INTRODUCTION · p. 2, RESULTS · p. 2

  • Results

    At low temperatures, PRIS finds ground states of arbitrary Ising models with high probability and probes their Gibbs distributions at varied temperatures.

    Source: Conclusion · p. 6

  • Takeaways & Limitations

    PRIS may help probe critical behavior despite critical slowing down by generating many samples per unit time, limited by active silicon photonics bandwidth.

    Source: Conclusion · p. 7

    A photonic realization still requires work to manage and characterize noise and to add homodyne detection for output amplitude and phase.

    Source: Conclusion · p. 7

Abstract

from arXiv · show

The inability of conventional electronic architectures to efficiently solve large combinatorial problems motivates the development of novel computational hardware. There has been much effort recently toward developing novel, application-specific hardware, across many different fields of engineering, such as integrated circuits, memristors, and photonics. However, unleashing the true potential of such novel architectures requires the development of featured algorithms which optimally exploit their fundamental properties. We here present the Photonic Recurrent Ising Sampler (PRIS), a heuristic method tailored for parallel architectures that allows for fast and efficient sampling from distributions of combinatorially hard Ising problems. Since the PRIS relies essentially on vector-to-fixed matrix multiplications, we suggest the implementation of the PRIS in photonic parallel networks, which realize these operations at an unprecedented speed. The PRIS provides sample solutions to the ground state of arbitrary Ising models, by converging in probability to their associated Gibbs distribution. By running the PRIS at various noise levels, we probe the critical behavior of universality classes and their critical exponents. In addition to the attractive features of photonic networks, the PRIS relies on intrinsic dynamic noise and eigenvalue dropout to find ground states more efficiently. Our work suggests speedups in heuristic methods via photonic implementations of the PRIS. We also hint at a broader class of (meta)heuristic algorithms derived from the PRIS, such as combined simulated annealing on the noise and eigenvalue dropout levels. Our algorithm can also be implemented in a competitive manner on fast parallel electronic hardware, such as FPGAs and ASICs.

INTRODUCTION

Heuristic methods face difficulty scaling to large combinatorial problems, while parallel photonic hardware offers a potential platform for fast computation. The paper introduces PRIS, a recurrent Ising sampler whose noisy updates and eigenvalue selection are analyzed for sampling and ground-state search.

INTRODUCTION

  • Heuristic methods struggle with large NP-hard optimization problems, motivating algorithms that can exploit parallel photonic hardware.The paper proposes PRIS as a fast heuristic based on iterative matrix multiplications.Source: INTRODUCTION · p. 2

RESULTS

  • PRIS encodes the spin state in parallel photonic signals, applies a linear transformation and a noisy nonlinear threshold, then feeds the updated state back.The proposed linear transform is derived from the eigendecomposition of the Ising coupling matrix.Source: RESULTS · p. 2, Figure 1 · p. 2

General theory of the PRIS dynamics

  • The recurrent update defines a transition probability; for symmetric couplings, it obeys detailed balance and supports an effective-Hamiltonian description.The expansion to the stated effective Hamiltonian is made in the large-noise approximation.Source: General theory of the PRIS dynamics · p. 3
  • Eigenvalue dropout level α controls how many coupling-matrix eigenvalues remain, with lower eigenvalues dropped when the real part of the square root is taken.The diagonal offset helps make the matrix under the square root positive definite for large positive α; the authors report that dropout can improve performance as α approaches zero.Source: General theory of the PRIS dynamics · p. 3
  • With Sqα(D) = D, PRIS matches the Ising free energy in the specified associative-memory setting, but uses synchronous rather than sequential spin updates.The authors note that synchronous updates are better suited to parallel photonic networks.Source: General theory of the PRIS dynamics · pp. 3–4
  • The general coupling regime is used for demonstrations, but additional considerations are needed to eliminate non-ergodic behavior.The passage notes both theoretical and simulation evidence for this concern.Source: General theory of the PRIS dynamics · p. 4

Finding the ground state of Ising models with the PRIS

The PRIS is evaluated as a heuristic ground-state search and as a sampler for critical behavior, with noise and eigenvalue dropout tuning its operation. The reported results connect these settings to ground-state search performance and show that PRIS samples recover properties of the two-dimensional Ising transition.

Finding the ground state of Ising models with the PRIS

  • On photonic parallel hardware, PRIS is projected to find ground states for graph orders N ∼100 in under a millisecond and remains roughly insensitive to graph density after noise optimization.The study also explores ground-state search in regimes around N ∼1,000 where exact solvers typically fail; static fabrication error remains relevant to fuller comparisons.Source: Finding the ground state of Ising… · p. 4

Influence of the noise and eigenvalue dropout levels

  • An optimal noise level balances declining ground-state probability at high noise against long autocorrelation times at low noise, minimizing iterations to find the ground state.The optimum can be approximated using ground-state probability and energy autocorrelation time.Source: Influence of the noise and eigenvalue… · p. 4
  • With α = 0, the tested problems had higher ground-state probability and lower autocorrelation times than with α > 0; some sparse models favored α < 0.Dropout perturbs the energy landscape and makes transition-matrix eigenvalues less uniform, while its autocorrelation-time reduction can counterbalance the perturbation; for the studied spin glass the ground state remained unchanged down to α = 0.Source: Influence of the noise and eigenvalue… · pp. 4–5

Detecting and characterizing critical behaviors with the PRIS

  • PRIS Binder-cumulant intersections estimated the 2D ferromagnetic Ising critical temperature as 2.241, within 1.3% of the theoretical 2.269.PRIS samples also exhibited the order-disorder transition, and estimated critical magnetization and susceptibility scaling parameters within 10% of theoretical values.Source: Detecting and characterizing critical… · p. 6
  • For critical sampling, parallel PRIS has estimated per-sample complexity O(N^(1+zPRIS/2)) versus O(N^(2+zMH/2)) for sequential MH, with an advantage when zPRIS < zMH + 2.The comparison concerns generating an uncorrelated sample; the PRIS estimate assumes a parallel architecture.Source: Detecting and characterizing critical… · p. 6, Figure 4 · p. 6

Conclusion

The conclusion presents PRIS as a heuristic sampler for Ising Gibbs distributions, with low-temperature ground-state search and implementations on photonic and electronic parallel hardware. It also highlights critical-behavior studies, possible dropout-based metaheuristics, and experimental requirements for a photonic realization.

Conclusion

  • An FPGA demonstration ran PRIS on N = 100 random spin glasses with 63 ns time steps, while the conclusion describes photonic clocks of ≲1 ns as a target.The paper frames the result as evidence that PRIS can use parallel architectures of different kinds, including hybrid opto-electronic implementations.Source: Conclusion · pp. 6–7
  • PRIS can probe Ising Gibbs distributions at finite temperatures and find ground states with high probability at low temperatures.The conclusion identifies matrix-to-vector accelerators including photonic networks, free-space optical processors, FPGAs, and ASICs.Source: Conclusion · pp. 6–7
  • The authors note that PRIS uses dynamic noise and suggest passivity-based computation and non-volatile phase-change materials as routes toward lower energy costs.The proposed energy-cost benefit is a possibility, not a demonstrated result in these passages.Source: Conclusion · p. 7
  • The conclusion proposes a broader photonic metaheuristic family, including simulated annealing over noise or eigenvalue dropout α.The authors suggest these variants could improve performance on larger graphs.Source: Conclusion · p. 7

ground states of Ising models

  • The reported ground-state benchmark uses random initial states and estimates the steps to reach a known ground state from run-length quantiles.Ground states were computed with BiqMac; the study also reports optimizing noise on a random sample graph for each graph-order and density pair.Source: ground states of Ising models · p. 7

HEURISTIC RECURRENT ALGORITHMS FOR PHOTONIC ISING MACHINES SUPPLEMENTARY INFORMATION

  • The supplementary overview covers theoretical analysis of PRIS dynamics, weight-matrix construction, numerical simulations, scaling, and comparisons with other heuristics.It describes the numerical work as estimating ground-state performance and sampling critical behavior.Source: HEURISTIC RECURRENT ALGORITHMS FOR… · p. 14

Definition of an Ising problem

  • The stated general Ising problem seeks a spin distribution minimizing a Hamiltonian defined by couplings and an external-field vector.The coupling matrix is symmetric; the work characterizes two subclasses, general antiferromagnets and spin glasses.Source: Definition of an Ising problem · p. 14

Mapping to the MAX-CUT problem

  • The paper maps general Ising problems to weighted MAX-CUT, allowing its ground-state solver to be used because the solutions coincide under K = −A.The MAX-CUT objective corresponds to dividing graph vertices into two spin groups to maximize edges between them.Source: Mapping to the MAX-CUT problem · p. 14

SUPPLEMENTARY NOTE 2: GENERAL THEORY OF PHOTONIC RECURRENT ISING SAMPLER

The algorithm is a synchronous, parallel recurrent network that updates neuron states from noisy matrix-multiplication potentials followed by a nonlinear threshold. Its temperature factor is calibrated by approximating transition probabilities with a sigmoid, with optimal scaling observed to depend linearly on noise scale.

Optical Parallel Recurrent MCMC algorithm

  • The framework generalizes noisy synchronous Hopfield dynamics to parallel recurrent sampling of hard Ising problems.Neurons and spins are used equivalently in the spin-glass setting.Source: Optical Parallel Recurrent MCMC… · p. 15
  • Each update applies a noisy matrix-multiplication potential to the neurons, then a nonlinear threshold to produce the next state.The algorithm starts from random neuron values and repeats the transformation; its variables are (C, fϕ, θ).Source: Optical Parallel Recurrent MCMC… · p. 15

Determination of the Temperature factor

  • Optimal sigmoid-rescaling factors vary linearly with noise scale, and Table I reports these factors and their corresponding errors across noise distributions.The analysis assumes symmetric noise; it notes an extension for infinite-variance distributions using a scaling parameter.Source: Determination of the Temperature factor · p. 16, Table I · p. 16

Transition probability

The transition probability is approximated by expanding around a sigmoid-based zero-order term, with the approximation error controlling higher-order contributions. The derivation treats symmetric noise and leaves nonsymmetric noise distributions for future extension.

Transition probability

  • The sigmoid approximation error serves as an expansion parameter, and higher-order terms are neglected when it is small.The text identifies ε ≪ s as the condition for the expansion.Source: Transition probability · p. 17, Expansion of the transition… · p. 17
  • The full transition probability is formed by multiplying independent single-spin flip probabilities, then expanded in powers of the approximation error ε.The zero-th-order term recovers the result of Ref.58, while the worst-case N-th-order term scales as ε^N.Source: Expansion of the transition… · p. 17
  • The derivation does not cover general nonsymmetric noise, though it proposes perturbation from the symmetric case or skew-normal parametrization as possible extensions.The skew-normal distribution is described as a natural extension of the analyzed symmetric toy model.Source: Expansion of the transition… · p. 17

Detailed balance condition

With symmetric couplings, the recurrent process satisfies detailed balance and has a stationary Gibbs distribution described by an effective Hamiltonian. The analysis develops small- and large-noise approximations, while mapping arbitrary Ising models in the large-noise regime remains incompletely characterized.

Detailed balance condition

  • The transition probability satisfies detailed balance when J is symmetric, enabling derivation of the effective Hamiltonian.The condition is expressed as equality of the cyclic sums of transition Hamiltonians.Source: Detailed balance condition · p. 18

Effective Gibbs distribution

Small noise approximation

Large noise approximation

  • The large-noise mapping of a general Ising Hamiltonian requires a symmetric real square root of K, and the general mapping study is left for future work.Complex-valued matrices are suggested as an interesting extension; the assumptions also include synchronous, Markovian dynamics with static connections.Source: Large noise approximation · p. 19
  • For the N = 16 random spin glass, the figure contrasts using all 16 eigenvalues with dropping negative eigenvalues to reduce the noise threshold of the large-noise expansion.The α = 0 case retains 7 of 16 eigenvalues.Source: Figure 1 · p. 20

Inequalities between Hamiltonians

Bounds on log cosh yield inequalities between Hamiltonians, which reverse direction when summed into partition-function inequalities. A norm-equivalence argument supplies another bound, and the Hamiltonians are illustrated for a sample random Ising problem.

Inequalities between Hamiltonians

Suggested algorithms

The PRIS uses a modified square root of the Ising coupling matrix to approximate the Gibbs distribution, with eigenvalue dropout and a diagonal offset shaping the search. The offset supports square-root construction, but diagonal dominance is not a necessary condition for positive definiteness.

Suggested algorithms

  • The suggested PRIS regime sets the recurrent loop C to a modified square root of K, with α controlling eigenvalue dropout and Δ a diagonal offset.An alternative is C = K, previously studied for associative-memory couplings.Source: Suggested algorithms · pp. 20–21

SUPPLEMENTARY NOTE 3: CONSTRUCTION OF THE WEIGHT MATRIX

  • Adding a diagonal offset makes the Ising matrix diagonally dominant and suitable for constructing a real symmetric square root used in the large-noise Hamiltonian.The construction adds Δ to K before taking its square root; the offset contributes only a spin-independent constant to the Hamiltonian.Source: SUPPLEMENTARY NOTE 3: CONSTRUCTION OF… · p. 22

Discussion on the diagonal offset

  • Diagonal dominance is sufficient but not necessary for positive definiteness, and the suggested offset is especially adapted to α > 0 and same-sign couplings.For large spin glasses, the authors typically use Δ_ii = |Σ_j≠i K_ij|; α can have a lower limit near −1/N for nonzero square-root matrices.Source: Discussion on the diagonal offset · pp. 22–23

Tuning the search dimensionality with eigenvalue dropout

Eigenvalue dropout changes the dimensionality of the PRIS search, with the heuristic suggesting that values around α = 0 can improve performance. A tunable-offset variant is reported to likely preserve the ground state when α ≥ 0.

Tuning the search dimensionality with eigenvalue dropout

  • Reducing α drops lower eigenvalues, and the heuristic suggests ground states may favor positive-eigenvalue eigenspaces, explaining observed performance improvements around α = 0.This is a heuristic explanation, not a guarantee: the hard constraint on the spin configuration must still hold.Source: Tuning the search dimensionality with… · p. 23

Modified algorithm with tunable offset

SUPPLEMENTARY NOTE 4: NUMERICAL SIMULATIONS

The numerical simulations examine PRIS sampling and critical-exponent estimation on two analytically solvable Ising models. They also describe limitations in computational scaling and temperature calibration, and compare critical-exponent fits with Metropolis-Hastings.

SUPPLEMENTARY NOTE 4: NUMERICAL SIMULATIONS

  • The simulations evaluate PRIS sampling of Gibbs distributions for square-lattice two-dimensional ferromagnetic and infinite-range Ising models, estimate critical exponents, and estimate the Binder cumulant as a function of temperature.Source: Evaluating sampling performance and… · pp. 25, 27
  • Adding a diagonal offset suppresses the odd dynamics seen with C = K, which can otherwise converge to the maximal-energy state in some runs.The figure caption says the offset version converges to the ground state after some time.Source: Evaluating sampling performance and… · p. 27, Figure 7 · p. 28
  • Taking the square root of the coupling matrix removes symmetry and sparsity reductions, while the proposed effective temperature T = kφ² does not match the theoretical 2D ferromagnetic critical temperature.The authors note that this makes critical-temperature estimation with PRIS problematic.Source: Evaluating sampling performance and… · p. 27
  • Critical-exponent fits compare PRIS with Metropolis-Hastings for susceptibility, energy autocorrelation time, magnetization, and magnetization autocorrelation time in both models.The supplied table caption identifies the table as a summary of fitted exponents but gives no exponent values.Source: Evaluating sampling performance and… · p. 28, Table II · p. 28, Figure 8 · p. 29

SUPPLEMENTARY NOTE 5: SCALING OF THE PRIS PERFORMANCE ON SEVERAL PHOTONIC ARCHITECTURES

The paper evaluates PRIS implementation on cascaded programmable Mach-Zehnder interferometers and discusses optical neural networks based on photoelectric multiplication. Phase-setting precision affects the simulated scaling, while the latter architecture is proposed for much larger graphs.

Cascaded arrays of programmable Mach-Zehnder interferometers (MZI)

  • A bV-bit voltage-setting accuracy corresponds approximately to (bV − 1)-bit phase accuracy, while bθ-bit phase accuracy corresponds to (bθ + 1)-bit voltage accuracy.The relation follows from the paper's quadratic approximation between phase and voltage.Source: Cascaded arrays of programmable… · p. 31
  • An 8-bit phase resolution may increase required PRIS steps by one to two orders of magnitude, whereas 16-bit resolution does not impact performance in the simulations.The impact varies with graph order and topology; the comparison concerns simulated ground-state search on cascaded MZI arrays.Source: Cascaded arrays of programmable… · pp. 31–32

Optical Neural Networks based on Photoelectric Multiplication

  • Photoelectric-multiplication networks are proposed for PRIS at N ∼ 10^3–10^6, with noise tunable through photons per multiply-accumulate and the nonlinear function executed electronically.For random spin glasses, the reported estimates are nmac ∼ 4 at N = 100 and ∼ 15 at N = 1,000; energy per multiplication is 6.2 ± 0.35 fJ and 1.9 ± 0.5 pJ, respectively.Source: Optical Neural Networks based on… · p. 32

Free space optical Neural Networks

The PRIS is compared with other heuristics and hardware options, including large-spin-glass benchmarks and the PRIS-A annealing variant. The section also outlines free-space photonic implementation and its practical limits.

Free space optical Neural Networks

  • Free-space optical neural networks28,40 using 3D-printed masks or reconfigurable SLMs are another option for PRIS at N ∼10^6 neurons; speed is limited by photodetector and modulation bandwidths, typically ∼10 THz and ∼1 kHz for SLMs, ⪆1 GHz.Source: Free space optical Neural Networks · p. 32

Comparison table of various heuristic algorithms and architectures

SUPPLEMENTARY NOTE 6: COMPARISON OF THE PRIS TO SEVERAL (META)HEURISTICS

  • Across 10 random spin glasses with N ∼1000, PRIS achieved mean ground-state energies similar to MH; BiqMac found only an approximate state after three hours.That approximate state was outperformed by MH, PRIS-A, and SA.Source: Benchmarking versus… · p. 34, Figure 11 · p. 34

PRIS-A: a proposed metaheuristic variation of the PRIS

  • PRIS-A anneals PRIS noise like simulated annealing and can also vary eigenvalue dropout α, but a comprehensive study of this algorithm class is beyond the work's scope.The authors also suggest annealing α to affect the dimensionality of the ground-state search.Source: PRIS-A: a proposed metaheuristic… · p. 35
  • For both SA and PRIS-A, λ > 0.98 yielded consistently low-energy ground states, with no particular improvement from increasing λ.The reported PRIS-A settings used φi = 50 and φf = 0.1.Source: PRIS-A: a proposed metaheuristic… · p. 35
  • MH found 0.53% lower ground-state energies than PRIS on average; PRIS-A improved on PRIS by 0.46%, while SA beat PRIS-A by 0.18%.The comparison used the settings described for Table IV.Source: PRIS-A: a proposed metaheuristic… · pp. 35–36, Table IV · p. 35

SUPPLEMENTARY NOTE 7: IMPLEMENTATION OF THE PRIS ON FPGA

The FPGA implementation uses a configurable memory-backed recurrent loop and a binary tree for matrix-vector multiplication. Its memory capacity constrains scaling, while tests indicate that limited-precision rounding can preserve algorithm performance.

Architecture

  • A ZCU104 FPGA implementation uses AXI to load matrix, noise, and threshold data into BRAM, run the recurrent loop, and return the final state.The design uses an address-mapped AXI interface and a finite-state-machine loop module.Source: Architecture · p. 37
  • Matrix-memory bandwidth is the scaling bottleneck: available memory is sufficient for N < 100, but becomes a major limitation as problem size grows.The board's BRAM and URAM provide a maximum combined data rate of 36,288 bits per clock cycle.Source: Architecture · p. 37
  • The binary-tree multiplier conditionally sums each matrix row using the state vector, with rows loaded onto tree leaves and propagated toward the root.The tree's propagation delay is determined by its height.Source: Architecture · p. 37, Figure 13 · p. 38
  • For N = 100, rounding matrix, noise, and threshold values after scaling by 32 or 64 reached essentially equivalent average energy to double precision.The comparison used 1,000 trials with different input matrices and 1,000 iterations per trial.Source: Architecture · p. 38, Figure 14 · p. 38

Complexity

The complexity analysis focuses on memory-bound matrix multiplication and tree propagation, then tests timing and correctness on an FPGA. The implementation matches the predicted per-step cycle count and reproduces Python-computed energies.

Complexity

  • Matrix multiplication dominates a step because the entire N × N matrix must be read; tree propagation adds O(log2 B) cycles, with B bounded by available data and row size.The analysis also accounts for limits on BRAM bandwidth and register capacity.Source: Complexity · p. 39

Experiment

  • Scaling below 64 showed deviation from the unrounded algorithm, while scaling above 1024 produced the same sequence of state vectors.The tests compare scaled, rounded inputs with the original floating-point run.Source: Experiment · p. 40
  • FPGA-simulated energies exactly matched Python-calculated energies, verifying the implementation's correctness for the tested run.The test loaded a scaled case onto the FPGA and output the state at each step for energy calculation.Source: Experiment · p. 40, Figure 16 · p. 40
  • The N = 100 FPGA implementation took approximately 63 ns per algorithm step, matching the predicted 19-cycle step at 300 MHz.The 19 cycles comprise 10 for matrix reading, 8 for binary-tree addition, and 1 for noise injection and thresholding.Source: Experiment · pp. 39–41

Related papers

Loading 1811.02705v3…