PaperlayerView exact revision on arXiv

Source-linked AI summary

Massively Parallel Probabilistic Computing with Sparse Ising Machines

Navid Anjum Aadit, Andrea Grimaldi, Mario Carpentieri, Luke Theogarajan, John M. Martinis, Giovanni Finocchio, Kerem Y. Camsari

arXiv:2110.02481v2cs.ETcond-mat.dis-nncs.DC

TL;DR

Computationally hard sampling and optimization are hindered by the serial updating required by standard Gibbs sampling. The paper introduces a sparsified, massively parallel Ising Machine using invertible-logic problem encodings and exact graph sparsification, reporting faster sampling, 32-bit factorization, and substantial SAT-solver speedups. The FPGA results are constrained by the number of p-bits fitting on one device, with larger implementations proposed through stronger FPGAs, analog CMOS, or nanodevices.

  • Problem

    Standard Gibbs sampling requires sequential updates, creating a fundamental serial bottleneck for probabilistic computing and hard optimization.

  • Method

    The sparse Ising Machine combines invertible-logic circuit encodings, graph coloring, exact sparsification, and massively parallel hardware for MCMC sampling and optimization.

  • Results

    The sIM outperforms CPU, GPU, TPU, and competition-winning SAT-solver baselines while reliably factoring semiprimes up to 32 bits.

  • Takeaways & Limitations

    Sparse problem representations enable near-ideal parallelism, with flips per second scaling linearly with the number of p-bits.

  • Takeaways & Limitations

    The FPGA experiments were limited by the number of probabilistic bits that fit on a single device.

Abstract

from arXiv · show

Inspired by the developments in quantum computing, building domain-specific classical hardware to solve computationally hard problems has received increasing attention. Here, by introducing systematic sparsification techniques, we demonstrate a massively parallel architecture: the sparse Ising Machine (sIM). Exploiting sparsity, sIM achieves ideal parallelism: its key figure of merit - flips per second - scales linearly with the number of probabilistic bits (p-bit) in the system. This makes sIM up to 6 orders of magnitude faster than a CPU implementing standard Gibbs sampling. Compared to optimized implementations in TPUs and GPUs, sIM delivers 5-18x speedup in sampling. In benchmark problems such as integer factorization, sIM can reliably factor semiprimes up to 32-bits, far larger than previous attempts from D-Wave and other probabilistic solvers. Strikingly, sIM beats competition-winning SAT solvers (by 4-700x in runtime to reach 95% accuracy) in solving 3SAT problems. Even when sampling is made inexact using faster clocks, sIM can find the correct ground state with further speedup. The problem encoding and sparsification techniques we introduce can be applied to other Ising Machines (classical and quantum) and the architecture we present can be used for scaling the demonstrated 5,000-10,000 p-bits to 1,000,000 or more through analog CMOS or nanodevices.

I. INTRODUCTION

The paper addresses the serial bottleneck of Gibbs sampling by combining invertible-logic sparsification with a massively parallel sparse Ising Machine architecture. The resulting sIM targets probabilistic sampling and hard combinatorial optimization, with reported gains over CPU, GPU, TPU, and SAT-solver baselines.

  • Sequential Gibbs updates serialize network evolution because parallel updates can cause repeated oscillations and prevent convergence to the Boltzmann distribution.
  • The sIM combines invertible logic, sparse graph representations, coloring, and parallel hardware to overcome sequential MCMC updating.The architecture implements chromatic Gibbs sampling, updating conditionally independent color blocks in parallel.
  • Exact sparsification preserves the optimization structure while reducing high-degree nodes, and the same ideas are proposed for other Ising Machines and hardware platforms.The paper notes that more sophisticated sampling algorithms and nanodevice implementations could provide additional speedup.
  • 4-700x runtime improvements over competition-winning SAT solvers were reported to reach 95% accuracy on Boolean satisfiability problems.
  • 5-18x measured speedup over optimized TPU and GPU implementations and up to 6 orders of magnitude over a standard CPU were reported for probabilistic sampling.
  • 32-bit semiprimes were reliably factored, exceeding reported results from D-Wave quantum annealers and similar probabilistic solvers.

III. COMPOSABILITY AND RECONFIGURABILITY

Invertible probabilistic gates compose Boolean functions into reconfigurable p-circuits for inverse problems. Graph coloring and sparsification then support hardware implementation across factorization and satisfiability instances.

  • Composability: Invertible logic uses probabilistic COPY/NOT, AND, OR, and full-adder gates to compose p-circuits representing Boolean functions.
  • Composability: Running invertible multiplier and Boolean circuits in reverse enables hardware solutions for integer factorization and satisfiability problems.
  • Reconfigurability: Reconfigurable circuits accommodate different factorization and SAT instances through bias clamping, input clamping, and multiplexing.
  • Reconfigurability: Reconfigurable sparsification requires approximately 2,000 spins for a 32-bit factorization graph versus approximately 10,000 spins for a Chimera target under minor graph embedding.
  • Reconfigurability: Maximum-Cut can be efficiently mapped to a Boolean satisfiability instance, adding another layer of reconfigurability.
  • Graph coloring: DSATUR colors sparse graphs so conditionally independent p-bits can be updated in parallel, while low-density graphs may use few colors despite irregular high-degree nodes.

8-bit Factorizer

The factorizer and satisfiability examples illustrate how graph coloring and exact sparsification support parallel sIM operation. Five-color scheduling enables parallel updates, while sparsification addresses high-degree hardware bottlenecks.

  • Factorizer architecture: The 8-bit factorizer graph is colored with 5 colors, and the architecture uses 5 parallel, equally phase-shifted clocks.
  • Graph sparsity: Exact sparsification reduces large neighbor counts at the expense of additional bits while preserving the optimization structure.
  • Factorizer architecture: Different color blocks receive different clocks so neighboring p-bits do not flip simultaneously during exact Gibbs sampling.
  • Factorizer architecture: The full network is updated in parallel within one clock period while preserving effectively sequential operation across color blocks.
  • Graph sparsity: Graph density ρ measures sparsity, with vertices representing p-bits and edges representing interconnections in the J matrix.
  • Graph sparsity: A 32-bit factorizer graph has 784 p-bits, density 1.03%, and requires 5 distinct colors.

B. Sparsification of problem graphs

The sIM removes high-degree bottlenecks through exact graph sparsification, trading additional p-bits for bounded neighbor counts and faster clocks. Its parallel FPGA implementation achieves linearly scaling flips per second and substantially lower factorization time than a serial CPU.

  • Sparsification motivation: High-degree p-bits create synapse delays that severely limit clock speed, motivating exact sparsification.The bottleneck arises from completing large weighted additions within a fraction of a clock period.
  • Sparsification method: COPY gates split one p-bit into two ferromagnetically coupled p-bits without changing the ground states of the optimization problem.The equivalence holds at the end of an annealing schedule with high β.
  • Factorization graph: 2128 p-bits and maximum degree k = 5 sparsify the 32-bit factorizer graph, reducing density to 0.2% while enabling fast clocks.The sparsified graph uses five colors and limits each p-bit to five neighbors.
  • 3SAT graph: 1935 p-bits and maximum degree k = 4 sparsify the 3SAT graph, reducing density to 0.2% and avoiding large adder delays.The original graph contains a 498-neighbor node where all clause outputs meet.
  • Performance scaling: 5.99×10^9 to 8.06×10^10 fps accompanies increasing problem size for the sIM, while CPU fps is limited to around 10^5 and declines with p-bit count.The sIM uses parallel p-bit updates with five phase-shifted 15 MHz clocks assigned by graph coloring.
  • Solution time: 2.4×10^6x improvement at 50-bit factorization results from sIM TTS99 of 1.84 ms versus 4408.93 s for the CPU.TTS99 measures time to reach 99% of the absolute ground energy; the sIM remains roughly constant from 14-bit to 50-bit problems.

B. GPU and TPU comparison

The sIM’s flips-per-second performance scales linearly with p-bit count and exceeds optimized GPU and TPU implementations in the reported comparisons. The architecture also supports projections to massively parallel nanodevice implementations, while exact factorization reaches 32-bit semiprimes but exhibits exponential time scaling.

  • 143.8 flips/ns at 4,793 p-bits, yielding 18.03x the single Nvidia Tesla C1060 GPU with multi-spin coding.
  • 11.17x and 12.65x more flips/ns than the Google Cloud TPU and reference Nvidia Tesla V100 GPU, respectively.
  • Linear flips-per-second scaling with p-bit count supports projected sIM architectures reaching 1 million flips per nanosecond with 10^6 MTJs, if connectivity remains sparse.
  • The sIM reliably factors random semiprimes up to 32 bits, reported as the largest factorization among compared hardware and probabilistic solvers.
  • Exact factorization time depends empirically exponentially on problem size, and the authors state this approach is impractical for cryptography.

VII. BOOLEAN SATISFIABILITY WITH INVERTIBLE LOGIC

The sIM applies invertible-logic encodings and sparse parallel updates to 3SAT, achieving ideal scaling and fastest approximate solutions, while exact solving remains difficult at larger sizes.

  • Encoding and architecture: The sIM encodes 3SAT instances as invertible probabilistic circuits and uses sparse graph structure for hardware acceleration.The architecture parallelizes updates across colored p-bit blocks.
  • Scaling performance: 1.44×10^11 fps was achieved on the largest instance with 4793 p-bits and 1065 clauses.Flips per second increased linearly with graph size, indicating ideal parallelism.
  • Exact solving: Exact solving was slower than professional SAT solvers, and the sIM did not reach the ground state beyond 2903 p-bits.The exact target required satisfying all clauses.
  • Approximate solving: 4 to 700 times faster approximate solutions were obtained than competition-winning SAT solvers at 95% clause satisfaction.The comparison included CDCL solvers and Yalsat95.
  • Inexact sampling: Moderate overclocking improved time to solution, but excessive message dropping caused a sharp decline in convergence performance.This behavior appeared in both factorization and 3SAT experiments.
  • Scaling outlook: The FPGA results are limited by the number of p-bits that fit on one device, while scaled CMOS or nanodevice implementations are projected to extend capacity.The proposed architecture is intended for larger implementations.

X. METHODS

The experiments generate factorization and 3SAT benchmarks, implement simulated annealing and timing control on an FPGA, and measure flips per second and time to solution against software baselines.

  • Benchmark generation: Random semiprimes from 14-bit to 50-bit were generated, with 10 different numbers per instance.The resulting factorization graphs were very sparse.
  • Benchmark generation: 3SAT instances were collected as CNF files from the UBC SATLIB library, with each clause containing exactly three variables.These graphs were also very sparse.
  • Annealing procedure: Simulated annealing gradually increases β, with updated J and h values sent from MATLAB to the FPGA over time.The FPGA performs the p-bit evolution using the updated parameters.
  • Measurement procedure: A programmable FPGA timer freezes all p-bits simultaneously before MATLAB reads their states through USB-JTAG.This supports controlled stopping and post-processing.
  • Performance measurement: Flips per second was measured from counted flips and elapsed time, with figure error bars based on 100 measurements.CPU measurements used MATLAB timing functions.
  • Baseline comparison: SAT solver runtimes were measured using online source code built on the same Linux machine, with CDCL solvers targeting 100% satisfaction.Yalsat reported time to its current best solution.

A. Characteristics of CMOS and nanodevice based p-bits

The paper describes p-bit implementations in CMOS, nanodevices, and FPGA hardware, including random generation, activation lookup, clocking, matrix computation, and binary-variable conversion.

  • CMOS and nanodevice p-bits: Digital CMOS implements a p-bit with an RNG, LUT, and comparator, while a 14 nm FinFET and stochastic MTJ provide a nanodevice alternative.Both implementations realize the p-bit update equation.
  • FPGA architecture: The FPGA implementation uses MATLAB communication, a clocking unit, and a MAC unit to operate the sparse Ising Machine.The design is implemented on a Xilinx Virtex UltraScale+ VCU118 board.
  • Parallel clocking: A 250 MHz LVDS system clock generates equally phase-shifted parallel clocks for p-bit blocks through the FPGA clock manager.The clocks drive LFSRs inside colored p-bit blocks.
  • MAC computation: The MAC unit uses 32-bit LFSRs as random-number generators and configures LUT width and comparison logic for p-bit updates.Several RNG families and bit widths were experimentally evaluated.
  • Variable representation: FPGA calculations use binary variables, converting bipolar representations so each output satisfies m_i ∈ {0, 1}.The activation LUT maps tanh to (1 + tanh)/2.

C. Basic logic gates for probabilistic computing

Basic probabilistic logic gates provide reusable building blocks for invertible circuits, while fusion and sparsification compose equivalent circuits with reduced node degree.

  • Basic gates: COPY, NOT, AND, and OR gates are represented with J and h matrices whose ground states match their Boolean truth tables.The energy plots correspond to the Boltzmann probabilities of the gate states.
  • Fusion: Fusion combines shared p-bit nodes and merges the component J and h matrices into larger composite-circuit matrices.An AND and OR composition changes [3 × 3] matrices into a [5 × 5] J matrix and a [5 × 1] h matrix.
  • Fusion: The fused circuit has ground states matching the corresponding Boolean truth table.The energy plot identifies the valid composite states.
  • Circuit equivalence: The sparse and fused circuits exhibit the same ground states, demonstrating equivalence of the two representations.The example splits m3 and introduces m′3 through a copy gate.
  • Sparsification: Sparsification splits one node into multiple nodes and inserts a copy gate to reduce vertex degree.The sparse circuit remains equivalent because the copy relation preserves the represented logic.

1. Mathematical equivalence of fused and sparsified graphs

As β increases, fused and sparsified circuits converge to the same ground-state probabilities because ferromagnetic coupling forces duplicated nodes to agree. This justifies composing and simplifying p-circuits through node fusion.

  • Ground-state equivalence: Fused and sparsified circuits therefore have identical final probabilities in the annealed limit for a given optimization problem.The comparison applies, for example, to the fused and sparsified versions of the same subcircuit.
  • Ground-state equivalence: As β →∞, the probability of disagreement between connected duplicate nodes approaches zero, independent of the coupling strength J_T.The disagreement probability is given as P(m_A ≠ m_B) = exp(−βJ_T).
  • Ground-state equivalence: Fusing nodes makes the coupling term a constant, so it cancels from the final Boltzmann probabilities.The result follows after substituting m_A = m_B and using m_A^2 = 1.
  • Fusion construction: The mathematical justification for fusion requires adding the columns and corresponding rows of fused nodes to preserve symmetry in J.This produces the composite circuit representation.
  • Fusion construction: A fused n-bit factorizer reduces the circuit to a generalized p-bit count built from AND-gate inputs, fused full-adder rows, and additional terms.The supplied expression is fact = 2m + (4m + 1) + (3m + 1)(m − 2) + 1.
  • Hardware trade-off: Although fusion keeps the state space smaller, excessive fan-out creates hardware issues and slows the sIM clock.These constraints motivate sparsification when some p-bits acquire too many neighbors.

3. Sparsified circuit

Sparsification limits each p-bit’s fan-out by adding copy-gate nodes, producing hardware-friendly circuits with shorter adder delays. Graph density decreases with problem size, while the sparsity choice trades p-bit count against flips per second.

  • Sparsification method: Sparsified p-circuits avoid fusing p-bits whenever fusion would exceed the prescribed maximum neighbor count k.The approach is demonstrated for integer factorization and 3SAT.
  • Sparsification method: Copy-gate trees split same-signal inputs so that each p-bit remains within the chosen fan-out bound k.For every k − 1 input bits, an additional p-bit is introduced recursively until one representative remains.
  • Circuit size: The unsparsified factorizer contains 3m^2 p-bits for AND gates and 5m(m − 1) p-bits for full adders.Sparsification adds further p-bits for copy gates.
  • Circuit size: The factorization experiments set the maximum neighbor count to k = 5.
  • Hardware impact: Limiting fan-out makes the sparsified circuit hardware-friendly by enabling faster clocks with smaller adder delays.
  • Performance trade-off: Under a fixed 100 W budget, the sparsest representation gives the fastest estimated flips per second, while denser graphs support further scaling at lower speed.The projections assume 20 µW per p-bit and consider up to one million p-bits.
  • Graph density: The maximum density of a graph with k neighbors and |V| nodes is approximately k/|V|.This follows from comparing k-regular and all-to-all edge counts.
  • Graph density: For both factorization and 3SAT, graph density decreases with problem size even without sparsification, using k = ∞.The passage concludes that these instances can be represented efficiently in sparse, scalable hardware.

E. Scalability analysis

The architecture scales linearly in p-bits and interconnections when maximum degree is fixed, but sparsity trades hardware resources and clock speed. In factorization, sIM embedding supports 32-bit instances with far fewer spins than the cited MGE comparison.

  • Resource scaling: O(N) scaling applies to both p-bit resources and interconnections when each p-bit has a fixed maximum neighbor count k.
  • Embedding comparison: A fixed 2128-spin sIM hardware can factor all integers up to 32 bits, whereas Chimera MGE requires approximately 10000 spins for the 32-bit factorizer.
  • Sparsity trade-off: The sparsest representations maximize flips per second but increase p-bit count, while clock frequency decreases linearly with k under the assumed adder-delay model.
  • Embedding comparison: MGE on the cited fixed topologies fails beyond 18 bits for King’s graphs and appears to fail for grid graphs.The comparison assumes approximately 50000-spin target hardware.
  • Embedding comparison: The sIM uses a k = 4 sparsified graph to factor 32-bit integers at all tested sizes.
  • Time to solution: An exponential fit over up to 5375 p-bits characterizes increasing factorization time beyond 32-bit instances under the current annealing schedule.The fitted parameters are t0 = 10^-3.39 s and τ = 122.13.
  • Time to solution: An exponential fit versus bit width up to 50 bits gives t0 = 10^-7.17 s and τ = 1.26.

G. Error models for inexact Gibbs sampling

The paper models overclocking-induced inexact Gibbs sampling with masks that alter neighbor updates. Simulations show both models remain close to the exact distribution at low error but diverge sharply beyond a threshold toward parallel-update behavior.

  • Model setup: The error analysis studies overclocking effects using two mask-based models applied to a five-p-bit full-adder circuit.The authors state that the conclusions and limits are generally applicable.
  • Mask models: The single-mask model uses +1 for functioning connections and −1 for failing connections.
  • Evaluation: Both models are tested across every error fraction from 0/20 to 20/20 using 400 random masks and 2 × 10^4 samples per mask.
  • Mask models: The double-mask model uses two complementary masks containing 1s and 0s for updated and non-updated values.
  • Distributional behavior: At low error fractions, both models show relatively small deviations from the exact Boltzmann distribution before worsening sharply near a threshold.
  • Distributional behavior: As the error fraction approaches one, the observed steady-state distribution approaches the parallel-update distribution.
Loading 2110.02481v2…