GPU-accelerated simulated annealing based on p-bits with real-world device-variability modeling

TL;DR

The paper tackles the impact of real-world device variability on p-bit–based simulated annealing (pSA) and demonstrates the need for a scalable, hardware-relevant simulator. It introduces a GPU-accelerated, open-source pSA framework that maps combinatorial problems to Ising energy $H(\sigma) = - sum_i h_i \sigma_i - sum_{i<j} J_{ij} \sigma_i \sigma_j$ with input $I_i(t) = I_0( h_i + sum_j J_{ij} \sigma_j(t) )$ and p-bits updated as $\sigma_i(t) = sgn( r_i(t) + tanh(I_i(t)) )$. It also models device variability via three parameters—timing $\nu_i$, intensity $\lambda_i$, and offset $\delta_i$—and uses CUDA to achieve about two orders of magnitude speedup over CPU on MAX-CUT (G-set) problems with 800–20,000 nodes, plus two pSA variations (TA pSA and SpSA). Key findings show that variability can both degrade and enhance performance, with TA pSA and SpSA offering robust, near-optimal solutions under realistic MTJ variability, providing a valuable toolkit for probabilistic computing research and hardware-oriented optimization.

Abstract

Probabilistic computing using probabilistic bits (p-bits) presents an efficient alternative to traditional CMOS logic for complex problem-solving, including simulated annealing and machine learning. Realizing p-bits with emerging devices such as magnetic tunnel junctions (MTJs) introduces device variability, which was expected to negatively impact computational performance. However, this study reveals an unexpected finding: device variability can not only degrade but also enhance algorithm performance, particularly by leveraging timing variability. This paper introduces a GPU-accelerated, open-source simulated annealing framework based on p-bits that models key device variability factors -- timing, intensity, and offset -- to reflect real-world device behavior. Through CUDA-based simulations, our approach achieves a two-order magnitude speedup over CPU implementations on the MAX-CUT benchmark with problem sizes ranging from 800 to 20,000 nodes. By providing a scalable and accessible tool, this framework aims to advance research in probabilistic computing, enabling optimization applications in diverse fields.

Figures (6)

• Figure 1: Simulated annealing using p-bits (pSA) operates based on the probabilistic nature of p-bits (top), as described by \ref{['eqn:pbits']}. A combinatorial optimization problem is mapped onto an Ising model, which corresponds to an energy function (Hamiltonian). In this model, each p-bit is biased by $h$ and interacts with other p-bits through weights $J$ (bottom left). pSA seeks to reduce the energy of the Ising model by altering the states of p-bits $\sigma_i$. When the global minimum energy $H_{min}$ is reached, the states $\sigma_i$ represent a solution to the combinatorial optimization problem (bottom right).
• Figure 2: Modeling of p-bits with device variability. The model introduces three new parameters to capture device variability: timing ($\nu_i$), intensity ($\lambda_i$), and offset ($\delta_i$). (a) These parameters accurately reflect the switching characteristics of each MTJ, enabling the model to account for variations in p-bit outputs due to changes in the MTJ dwell time. (b) The timing parameter $\nu_i$ captures shifts in the output signal in response to input signal changes. (c) The intensity parameter $\lambda_i$ captures the variability in the steepness of the $\tanh(\lambda_i I_i)$ function, indicating intensity to input changes. (d) The offset parameter $\delta_i$ captures shifts in the input threshold, which affect the transition points of the output curve. These parameters combined allow the model to reflect probabilistic behavior variations in p-bits according to the properties of each MTJ.
• Figure 3: Performance comparison of pSA, TApSA, and SpSA algorithms on CPU and GPU across different problem sizes without the device variability. (a) Annealing time as a function of the number of nodes: GPU implementations show significantly faster annealing times compared to their CPU counterparts, while pSA (GPU) exhibits longer times for larger problem sizes. (b) Normalized mean cut value as a function of the number of nodes: The SpSA and TApSA algorithms maintain high normalized mean cut values across all problem sizes, indicating their robustness in achieving optimal solutions. (c) Mean cut value versus annealing time for different algorithms in G1: As the annealing time increases, TApSA and SpSA converge to higher mean cut values, particularly in the GPU implementation. pSA struggles to find higher cut values efficiently on both CPU and GPU. (d) Extended analysis of mean cut value versus annealing time in G81: The GPU-based implementations of TApSA and SpSA consistently achieve higher mean cut values faster than their CPU implementations. In contrast, pSA on the GPU lags in performance relative to the other algorithms.
• Figure 4: Impact of variability parameters on mean cut value for different simulated annealing algorithms for G1. Each subplot shows the sensitivity analysis of three algorithms (pSA, TApSA, and SpSA) under variations in key variability parameters: $\sigma_{\lambda}$, $\sigma_{\delta}$, and $\sigma_{\nu}$. (a), (b), and (c): Effects of $\sigma_{\lambda}$ under different combinations of $\sigma_{\delta}$ and $\sigma_{\nu}$. When both $\sigma_{\delta} = 0$ and $\sigma_{\nu} = 0$, the mean cut values for pSA improve steadily with increasing $\sigma_{\lambda}$, while TApSA and SpSA maintain consistently high values (a). As $\sigma_{\delta}$ increases to 0.5 (b), the effect on pSA’s mean cut values diminishes slightly. The trends remain stable when $\sigma_{\nu}$ increases to 0.5 (c), showing resilience in TApSA and SpSA. Notably, in case (c), pSA also achieves a high mean cut value, indicating its robustness under these conditions. (d), (e), and (f): Effects of $\sigma_{\delta}$ under different combinations of $\sigma_{\lambda}$ and $\sigma_{\nu}$. With $\sigma_{\lambda} = 0$ and $\sigma_{\nu} = 0$ (d), the pSA shows a gradual improvement in mean cut values with increasing $\sigma_{\delta}$, while TApSA and SpSA maintain stable high performance. When $\sigma_{\lambda}$ increases to 0.5 (e), the patterns remain similar. As $\sigma_{\nu}$ increases to 0.5 (f), the stability in SpSA and TApSA continues. (g), (h), and (i): Effects of $\sigma_{\nu}$ under different combinations of $\sigma_{\lambda}$ and $\sigma_{\delta}$. For $\sigma_{\lambda} = 0$ and $\sigma_{\delta} = 0$ (g), the impact on pSA is evident, with an initial increase in mean cut values. As $\sigma_{\lambda}$ increases to 0.5 (h), the stability of TApSA and SpSA remains unaffected, while pSA shows a slight decrease in mean cut values. When $\sigma_{\delta}$ increases to 0.5 (i), all algorithms demonstrate high stability, with SpSA consistently performing the best.
• Figure 5: Impact of variability in $\sigma_{\lambda}$, $\sigma_{\delta}$, and $\sigma_{\nu}$ on the normalized mean cut value for different node sizes. Each subplot shows how the normalized mean cut value changes with increasing number of nodes under various levels of variability for pSA, TApSA, and SpSA. (a): $\sigma_{\lambda} = 0.5$, $\sigma_{\delta} = 0$, and $\sigma_{\nu} = 0$: SpSA and TApSA maintain high normalized mean cut values across all node sizes, while pSA shows lower cut values with a significant drop at higher node counts. (b): $\sigma_{\lambda} = 1$, $\sigma_{\delta} = 0$, and $\sigma_{\nu} = 0$: As $\sigma_{\lambda}$ increases to 1, SpSA continues to perform well, but TApSA shows a small decrease in performance, and pSA remains low. (c): $\sigma_{\lambda} = 0$, $\sigma_{\delta} = 0.5$, and $\sigma_{\nu} = 0$: Variability in $\sigma_{\delta}$ affects pSA more strongly, reducing its cut values, while TApSA and SpSA remain largely unaffected. (d): $\sigma_{\lambda} = 0$, $\sigma_{\delta} = 1$, and $\sigma_{\nu} = 0$: Further increase in $\sigma_{\delta}$ results in little to no effect on SpSA and TApSA, while pSA still performs poorly. (e): $\sigma_{\lambda} = 0$, $\sigma_{\delta} = 0$, and $\sigma_{\nu} = 0.5$: Introducing variability in $\sigma_{\nu}$ shows resilience in both SpSA and TApSA across all node sizes, but pSA remains significantly affected. (f): $\sigma_{\lambda} = 0$, $\sigma_{\delta} = 0$, and $\sigma_{\nu} = 1$: Even under high variability in $\sigma_{\nu}$, SpSA and TApSA maintain high normalized cut values, whereas pSA remains sensitive and continues to increase the normalized mean cut values.