CMOS plus stochastic nanomagnets enabling heterogeneous computers for probabilistic inference and learning

Extending Moore’s law by augmenting complementary-metal-oxide semiconductor (CMOS) transistors with emerging nanotechnologies (X) has become increasingly important. One important class of problems involve sampling-based Monte Carlo algorithms used in probabilistic machine learning, optimization, and quantum simulation. Here, we combine stochastic magnetic tunnel junction (sMTJ)-based probabilistic bits (p-bits) with Field Programmable Gate Arrays (FPGA) to create an energy-efficient CMOS + X (X = sMTJ) prototype. This setup shows how asynchronously driven CMOS circuits controlled by sMTJs can perform probabilistic inference and learning by leveraging the algorithmic update-order-invariance of Gibbs sampling. We show how the stochasticity of sMTJs can augment low-quality random number generators (RNG). Detailed transistor-level comparisons reveal that sMTJ-based p-bits can replace up to 10,000 CMOS transistors while dissipating two orders of magnitude less energy. Integrated versions of our approach can advance probabilistic computing involving deep Boltzmann machines and other energy-based learning algorithms with extremely high throughput and energy efficiency.

With the slowing down of Moore's Law [1], there has been a growing interest in domain-specific hardware and architectures to address emerging computational challenges and energy-efficiency, particularly borne out of machine learning and AI applications.One promising approach is the co-integration of traditional complementary metaloxide semiconductor (CMOS) technology with emerging nanotechnologies (X), resulting in CMOS + X architectures.The primary objective of this approach is to augment existing CMOS technology with novel functionalities, by enabling the development of physics-inspired hardware systems that realize energy-efficiency, massive parallelism, and asynchronous dynamics, and apply them to a wide range of problems across various domains.
Being named one of the top 10 algorithms of the 20 th century [2], Monte Carlo methods have been one of the most effective approaches in computing to solve computationally hard problems in a wide range of applications, from probabilistic machine learning, optimization to quantum simulation.Probabilistic computing with p-bits [3] has emerged as a powerful platform for executing these Monte Carlo algorithms in massively parallel [4,5] and energyefficient architectures.p-bits have been shown to be applicable to a large domain of computational problems from combinatorial optimization to probabilistic machine learning and quantum simulation [6][7][8].
Several p-bit implementations that use the inherent stochasticity in different materials and devices have been proposed, based on diffusive memristors [9], resistive RAM [10], perovskite nickelates [11], ferroelectric transistors [12], single photon avalanche diodes [13], optical parametric oscillators [14] and others.Among alternatives sMTJs built out of low-barrier nanomagnets have demonstrated significant potential due to their ability to amplify noise, converting millivolts of fluctuations to hundreds of millivolts over resistive networks [15], unlike alternative approaches with amplifiers [16].Another advantage of sMTJ-based pbits is the continuous generation of truly random bitstreams without the need to be reset in synchronous pulsebased designs [17,18].
The possibility of designing energy-efficient p-bits using low-barrier nanomagnets has stimulated renewed interest in material and device research with several exciting demonstrations from nanosecond fluctuations [19][20][21] to better theoretical understanding of nanomagnet physics [22][23][24][25] and novel magnetic tunnel  5) Ta (5) Asynchronous Clock The proposed sMTJ-based p-bit circuit with two branches whose outputs are provided to an operational amplifier.R ave is the average resistance of R P and R AP of the sMTJ.5 sMTJ-based p-bits provide tunable, truly random and asynchronous clocks to a digital field programmable gate array (FPGA).(c) Digital p-bits in the FPGA use lookup tables (LUT), comparators, synaptic weights, and pseudorandom number generators (PRNG).The clocks of the PRNG are driven by the truly random asynchronous outputs coming from the analog p-bits.(d) Pictorial representation of perpendicular sMTJ.(e) Image of a single p-bit circuit.(f) Image of the FPGA.The asynchronous clocks are input through the peripheral module (PMOD) pins.(g) Typical output of p-bits #1 to 5 using 5 sMTJs obtained from the p-bit circuit (see Supplementary Section III), showing variations in fluctuations.(h) Experimentally measured ⟨V OUT ⟩ the time average (over a period of 3 minutes) of the p-bit circuit output, as a function of DC input voltage V IN .The yellow squares are experimental data, and the blue dashed line is a fit of the form junction designs [26,27].Despite promising progress with hardware prototypes [28][29][30][31][32], large-scale probabilistic computing using stochastic nanodevices remains elusive.As we will establish in this paper, designing purely CMOS-based high-performance probabilistic computers suited to sampling and optimization problems is prohibitive beyond a certain scale (>1M p-bits) due to the large area and energy costs of pseudorandom number generators.As such, any large-scale integration of probabilistic computing will involve strong integration with CMOS technology in the form of CMOS+X architectures.Given the unavoidable device-to-device variability, the interplay between continuously fluctuating stochastic nanodevices (e.g., sMTJs) with deterministic CMOS circuits and possible applications of such hybrid circuits remain unclear.
In this paper, we first introduce the notion of a heterogeneous CMOS+sMTJ system where the asynchronous dynamics of sMTJs control digital circuits in a standard CMOS Field Programmable Gate Array (FPGA).We view the FPGA as a "drop-in replacement" for eventual integrated circuits where sMTJs could be situated on top of CMOS.Unlike earlier implementations where sMTJs were primarily used to implement neurons and CMOS or analog components circuits for synapses [28,29], we design hybrid circuits where sMTJ-based p-bits control a large number of digital circuits residing in the FPGA without dividing the system into neurons (sMTJ) and synapses (CMOS).We show how the true randomness injected into deterministic CMOS circuits augment low-quality random number generators based on linear feedback shift registers (LFSR).This result represents an example of how sMTJs could be used to reduce footprint and energy consumption in the CMOS underlayer.In this work, we present a small example of a CMOS + sMTJ system, however, similar systems can be scaled up to much bigger densities, leveraging the proven manufacturability of magnetic memory at gigabit densities.Our results will help lay the groundwork for larger implementations in the presence of unavoidable device-to-device variations.We also focus beyond the common use case of combinatorial optimization of similar physical computers [33], considering probabilistic inference and learning in deep energy-based models.
Specifically, we use our system to train 3-hidden 1visible layer deep and unrestricted Boltzmann machines that entirely rely on the asynchronous dynamics of the stochastic MTJs.Second, we evaluate the quality of randomness directly at the application level through probabilistic inference and deep Boltzmann learning.This approach contrasts with the majority of related work, which typically conducts statistical tests at the single device level to evaluate the quality of randomness [21,[34][35][36][37][38] (see Supplementary Sections VIII, XI, and XII for more randomness experiments).As an important new result, we find that the quality of randomness matters in machine learning tasks as opposed to optimization tasks that have been explored previously.And finally, we conduct a comprehensive benchmark using an experimentally calibrated 7-nm CMOS PDK and find that when the quality of randomness is accounted for, the sMTJbased p-bits are about 4 orders of magnitude smaller in area and they dissipate 2 orders of magnitude less energy, compared to CMOS p-bits.We envision that large-scale CMOS+X p-computers (≫ 10 5 ) can be a reality in scaled up versions of the CMOS + sMTJ type computers we discuss in this work.

Constructing the heterogeneous p-computer
FIG. 1 shows a broad overview of our sMTJ-FPGA setup along with device and circuit characterization of sMTJ p-bits.Unlike earlier p-bit demonstrations with sMTJs as standalone stochastic binary neurons, in this work, we use sMTJ-based p-bits to generate asynchronous and truly random clock sources to drive digital p-bits in the FPGA (FIG.1a,b,c).
The conductance of the sMTJ depends on the relative angle θ between the free and the fixed layers, G MTJ ∝ [1 + P 2 cos(θ)], where P is the interfacial spin polarization.When the free layer is made out of a low barrier nanomagnet θ becomes a random variable in the presence of thermal noise, causing conductance fluctuations between the parallel (P) and the antiparallel (AP) states (FIG.1d).
The five sMTJs used in the experiment are designed with a diameter of 50 nm and have a relaxation time of about 1 to 20 ms, with energy barriers of ≈14-17 k B T , assuming an attempt time of 1 ns [40] (see Supplementary Section II).In order to convert these conductance fluctuations into voltages, we design a new p-bit circuit (FIG.1b,e).This circuit creates a voltage comparison between two branches controlled by two transistors, fed to an operational amplifier.As we discuss in Supplementary Section III, the main difference of this circuit compared to the earlier 3 transistor/1MTJ design used in earlier demonstrations [28,29] is in its ability to provide a larger stochastic window to tune the p-bit (FIG.1h) with more variation tolerance (see Supplementary Section IV).
FIG. 1f,g show how the asynchronous clocks obtained from p-bits with 50/50 fluctuations are fed to the FPGA.Inside the FPGA, we design a digital probabilistic computer where a pbit includes a lookup table (LUT) for the hyperbolic tangent function, a pseudorandom number generator (PRNG) and a digital comparator (see Supplementary Section V).
The crucial link between analog p-bits and the digital FPGA is established through the clock of the PRNG used in the FPGA, where a multitude of digital p-bits can be asynchronously driven by analog p-bits.As we discuss in Sections 3-4, depending on the quality of the chosen PRNG, the injection of additional entropy through the clocks has a considerable impact on inference and learning tasks.The potential for enhancing low-quality PRNGs using compact and scalable nanotechnologies, such as sMTJs, which can be integrated as a BEOL (Back-End-Of-Line) process on top of the CMOS logic, holds significant promise for future CMOS + sMTJ architectures.

Probabilistic inference with heterogeneous p-computers
In the p-bit formulation, we define probabilistic inference as generating samples from a specified distribution which is the Gibbs-Boltzmann distribution for a given network (see Supplementary Section I for details).This is a computationally hard problem [41], and is at the heart of many important applications involving Bayesian inference [42], training probabilistic models in machine learning [43], statistical physics [44] and many others [45].Due to the broad applicability of probabilistic inference, improving key figures of merit such as probabilistic flips per second (sampling throughput) and energy-delay product for this task are extremely important.
To demonstrate this idea, we evaluate probabilistic inference on a probabilistic version of the full adder (FA) [39] as shown in FIG.2a.The truth table of the FA is given in FIG.2b.The FA performs 1-bit binary addition and it has three inputs (A, B, Carry in=C in ) and two outputs (Sum=S, and Carry out=C out ).The probabilistic FA can be described in a 5 p-bit, fully-connected network (FIG.2a).When the network samples from its equilibrium, it samples states corresponding to the truth table, according to the Boltzmann distribution.We demonstrate probabilistic sampling on the probabilistic FA using the digital p-bits with standalone Linear Feedback Shift Registers LFSRs (only using the FPGA), sMTJ-clocked LFSRs (using sMTJ-based p-bits and the FPGA), and standalone Xoshiro RNGs (only using the FPGA).Our main goal is to compare the quality of randomness obtained by inexpensive but low-quality PRNGs such as LFSRs [46] with sMTJ-enhanced LFSRs and high-quality but expensive PRNGs such as Xoshiro [47] (see Supplementary Section VI).
FIG. 2c shows the comparison of these three different solvers where we measure the Kullback-Leibler (KL) divergence [48] between the cumulative distribution based on the number of sweeps and the ideal Boltzmann distribution of the FA: where P exp is the probability obtained from the experiment (cumulatively measured) and P ideal is the probability obtained from the Boltzmann distribution.For LFSR (red line), the KL divergence saturates when the number of sweeps exceeds N = 10 4 , while for sMTJ-clocked LFSR (blue line) and Xoshiro (green line), the KL divergence decreases with increasing the number of sweeps.The persistent bias of the LFSR is also visible in the partial histogram of probabilities measured at N = 10 6 sweeps as shown in FIG.2d (see Supplementary Section VII for the full histograms).It is important to note here in our present context where sMTJs are limited to a handful of devices, we use sMTJbased p-bits used to drive low-quality LFSRs, observing how they perform similar to high-quality PRNGs.In integrated implementations however, sMTJ-based p-bits can be directly used as p-bits themselves, without any supporting PRNG (see Supplementary Section XVI for details on projections of integrated implementations).The mechanism of how the sMTJ-clocked LFSRs produce random numbers is interesting: even though the next bit in an LFSR is always perfectly determined, the randomness in the arrival times of clocks from the sMTJs makes their output unpredictable.Over the course of the full network's evolution, each LFSR produces an unpredictable bitstream, functioning as truly random bits.
The observed bias of the LFSR can be due to several reasons: first, the LFSRs generally provide low-quality random numbers and do not pass all the tests in the NIST statistical test suite [49] (see Supplementary Section XII).Second, we take whole words of random bits from the LFSR to generate large random integers.This is a known danger when using LFSRs [50,51], which can be mitigated by the use of phase shifters that scramble the parallelly obtained bits to reduce their correlation [52].But such measures increase the complexity of PRNG designs further limiting the scalable implementation of digital p-computers (see Supplementary Section XI for detailed experimental analysis of LFSR bias).
The quality of randomness in Monte Carlo sampling is a rich and well-studied subject (see, for example, [53][54][55]).The main point we stress in this work is that even compact and inexpensive simple PRNGs can perform as well as sophisticated, high-quality RNGs when augmented by truly random nanodevices such as sMTJs.

Boltzmann Learning with heterogeneous p-computers
We now show how to train deep Boltzmann machines (DBM) with our heterogeneous sMTJ + FPGA-based computer.Unlike probabilistic inference, in this setting, the weights of the network are unknown and the purpose of the training process is to obtain desired weights for a given truth table, such as the full adder (see Supplementary Section IX for an example of arbitrary distribution generation using the same learning algorithm).We consider this demonstration as a proof-of-concept for eventual larger-scale implementations (FIG.3a,b).Similar to probabilistic inference, we compare the performance of three solvers: LFSR-based, Xoshiro-based and sMTJ+LFSR-based RNGs.We choose a 32-node Chimera lattice [56] to train a probabilistic full adder with 5 visible nodes and 27 hidden nodes in a 3-layer DBM (see FIG. 3b top panel).Note that this deep network is significantly harder to train than training fully-visible networks whose data correlations are known a priori [29], necessitating positive and negative phase computations (see Supplementary Section VII and Algorithm 1 for details on the learning algorithm and implementation).
FIG. 3c,d show the KL divergence and the probability distribution of the full adder Boltzmann machines based on the fully digital LFSR/Xoshiro and the heterogeneous sMTJclocked LFSR RNGs.The KL divergence in the learning experiment is performed like this: after each epoch during training, we save the weights in the classical computer and perform probabilistic inference to measure the KL distance between the learned and ideal distributions.The sMTJclocked LFSR and the Xoshiro based Boltzmann machines produce probability distributions that eventually closely approximate the Boltzmann distribution of the full adder.On the other hand, the fully digital LFSR based Boltzmann machine produces the incorrect states [A B C in S C out ] = 2 and 29 with a significantly higher probability than the correct peaks, and grossly underestimates the probabilities of states 0, 6, 18, and 25 (see FIG. Supplementary Figure 4 for full histograms that are avoided here for clarity).As in the inference experiment (FIG.2a), the KL divergence of the LFSR saturates and never improves beyond a point.The increase in the KL divergence for Xoshiro and sMTJ-clocked LFSR towards the end is related to hyperparameter selection and unrelated to RNG quality [57].For this reason, we select the weights at epoch=400 for testing to produce the histogram in FIG.3d.
In line with our previous results, the learning experiments confirm the inferior quality of LFSR-based PRNGs, particularly for learning tasks (see Supplementary Section X for an MNIST training comparisons between p-bits based on Xoshiro and LFSR).While LFSRs can produce correct peaks with some bias in optimization problems, they fail to learn appropriate weights for sampling and learning, rendering them unsuitable for these applications.In addition to these results, statistical tests on the NIST test suite corroborate our findings that sMTJ-clocked LFSRs and high-quality PRNGs such as Xoshiro outperform the pure LFSR-based p-bits (see Supplementary Section XII).
Our learning result demonstrates how asynchronously interacting p-bits can creatively combine with existing CMOS technology.Scaled and integrated implementations of this concept could lead to a resurgence in training powerful deep Boltzmann machines [58].

Energy and transistor count comparisons
Given our prior results stressing how the quality of randomness can play a critical role in probabilistic inference and learning, it is beneficial to perform precise, quantitative comparisons with the various digital PRNGs we built in hardware FPGAs with sMTJ-based p-bits [15].Note that for this comparison, we do not consider augmented CMOS p-bits, but directly compare sMTJ-based mixed signal pbits with their digital counterparts (see Supplementary Section XVI for details on projections of integrated implementations using sMTJ-based mixed signal p-bits).Moreover, instead of benchmarking the voltage comparator based p-bit circuit shown in FIG. 1 or other types of spinorbit torque based p-bits [3,59], we benchmark the 3T/1MTJ based p-bit first reported in [15].The reason for this choice is that this design allows the use of fast in-plane sMTJs whose fluctuations can be as fast as micro to nanoseconds.We also note that the table-top components we use in this work are not optimized but used for convenience.
For the purpose of benchmarking and characterization, we synthesize circuits for LFSR and Xoshiro PRNGs and these PRNG-based p-bits using the ASAP 7nm Predictive process development kit (PDK) that uses SPICE-compatible FinFET device models [60].Our synthesis flow, explained in detail in Supplementary Section XIII, starts from hardware description level (HDL) coding of these PRNGs and leads to transistor-level circuits using the experimentally benchmarked ASAP 7nm PDK.As such, the analysis we perform here offers a high degree of precision in terms of transistor counts and quantitative energy consumption.FIG.4a shows the transistor count for p-bits using 32bit PRNGs.Three pieces make up a digital p-bit: PRNG, LUT (for the activation function) and a digital comparator (typically small).To understand how each piece contributes to the transistor count, we separate the PRNG from the LUT contributions in FIG.4a.
First, we reproduce earlier results reported in Ref. [28], corresponding to the benchmarking of the design reported in [15] and find that a 32-bit LFSR requires 1122 transistors which is very close to the custom-designed 32-bit LFSR with 1194 transistors in Ref. [28].However, we find that the addition of a LUT, ignored in [28], adds significantly more transistors.Even though the inputs to the p-bit are 10-bits (s [6][3]), the saturating behavior of the tanh activation allows reductions in LUT size.In our design, the LUT stores 2 8 words of 32-bit length that are compared to the 32-bit PRNG.Under this precision, the LUT increases the transistor count to 5150, and more would be needed for finer representations.Note that the compact sMTJ-based p-bit design proposed in [15] uses 3 transistors plus an sMTJ which we estimate as Energy per Random bit (fJ) Fig. 4: Transistor counts and energy consumption for p-bit and RNG implementations The digital p-bits and PRNGs are synthesized by the ASAP7 PDK and simulated in HSPICE in transistor-level simulations.The PRNGs are 32-bits long and LUTs store 2 8 words that are 32-bits long to be compared with 32-bit RNGs.The sMTJ-based p-bit result is repeated from [28].To activate the LUT, a periodic input signal with low inputs to the p-bit has been used.See the text and Supplementary Section XV for details on the energy calculation.
having an area of 4 transistors, following Ref.[28].In this case, there is no explicit need for a LUT or a PRNG.Additionally, the results presented in FIG. 2 and 3 indicate that to match the performance of the sMTJ-based p-bits, more sophisticated PRNGs like Xoshiro must be used.In this case, merely the PRNG cost of a 32-bit Xoshiro is 7516 transistors.The LUT costs are the same as LFSR-based p-bits which is about ≈ 4029 transistors.
Collectively, these results indicate that to truly replicate the performance of an sMTJ-based p-bit, the actual transistor cost of a digital design is about 11,000 transistors which is an order of magnitude worse than the conservative estimation performed in Ref. [28].
In FIG.4b we show the energy costs of these differences.We focus on the energy required to produce one random bit.Once again, our synthesis flow, followed by ASAP7 based HSPICE simulations reproduces the results presented in Ref. [28].We estimate a 23 fJ energy per random bit from the LFSR based PRNG where this number was reported to be 20 fJ in Ref. [28].
Similar to the transistor count analysis, we consider the effect of the LUT on the energy consumption which was absent in [28].We first observe that if the LUT is not active, i.e., if the input I i to the p-bit is not changing, the LUT does not change the energy per random bit very much.In a real p-circuit computation, however, I i would be continuously changing activating the LUT repeatedly.To simulate these working conditions, we create a variable I i pulse that wanders around the stochastic window of the pbit by changing the least significant bits of the input (see Supplementary Section XV).We choose a 1 GHz frequency for this pulse mimicking an sMTJ with a lifetime of 1 ns.We observe that in this case the total energy to create a random bit on average increases by a factor of 6× for the LFSR, reaching 145 fJ per bit.
For the practically more relevant Xoshiro, the average consumption per random bit reaches around 293 fJ.Once again, we conclude that the 20 fJ per random bit, reported in Ref. [28] underestimates the costs of RNG generation by about an order of magnitude when the RNG quality and other peripheries such as LUTs are carefully taken into account.In this paper, we do not reproduce the energy estimation of the sMTJ-based p-bit but report the estimate in Ref. [28] which assumes an sMTJ-based p-bit with ≈ nanosecond fluctuations.
Our benchmarking results highlight the true expense of high-quality, digital p-bits in silicon implementations.Given that functionally interesting and sophisticated p-circuits require above 10000 to 50000 p-bits [5], using a 32-bit Xoshiro-based p-bit in a digital design would consume up to 0.1 to 0.5 Billion transistors, just for the p-bits.In addition, the limitation of not being able to parallelize or fit more random numbers in hardware would limit the throughput [61] and the probabilistic flips per second, a key metric measuring the effective sampling speed of a probabilistic computer (see for example, [62][63][64]).As discussed in detail in Supplementary Section XVI, near-term projections with N = 10 4 p-bits using sMTJs with in-plane magnetic anisotropy (IMA) (τ ≈ 1ns [19]) can reach ≈ 10 4 flips/ns in sampling throughput.These results clearly indicate that a digital solution beyond 10000 to 50000 p-bits, as required by large-scale optimization, probabilistic machine learning, and optimization tasks, will remain prohibitive.To solve these traditionally expensive but practically useful problems, the heterogeneous integration of sMTJs holds great promise both in terms of scalability and energy-efficiency.

DISCUSSIONS
This work demonstrates the first hardware demonstration of a heterogeneous computer combining versatile FPGAs with stochastic MTJs for probabilistic inference and deep Boltzmann learning.We introduce a new variation tolerant p-bit circuit that is used to create an asynchronous clock domain, driving digital p-bits in the FPGA.In the process, the CMOS + sMTJ computer shows how commonly used and inexpensive PRNGs can be augmented by magnetic nanodevices to perform as well as high-quality PRNGs (without the resource overhead), both in probabilistic inference and learning experiments.Our CMOS + sMTJ computer also shows the first demonstration of training a deep Boltzmann network in a 32-node Chimera topology, leveraging the asynchronous dynamics of sMTJs.Careful comparisons with existing digital circuits show the true potential of integrated sMTJs which can be scaled up to million p-bit densities far beyond the capabilities of present day CMOS technology (see Supplementary Section XVI for detailed benchmarking and a p-computing roadmap).

METHODS sMTJ fabrication and circuit parameters
We employ a conventional fixed and free layer sMTJ, both having perpendicular magnetic anisotropy.The reference layer thickness is 1 nm (CoFeB) while the free layer is 1.8 nm (CoFeB), deliberately made thicker to reduce its energy barrier [28,35].The stack structure of the sMTJs we use is, starting from the substrate side, Ta(5)/ Pt(5)/ [Co(0.4)/Pt(0.4)]6 / Co(0.4)/Ru(0.4)/ [Co(0.4)/Pt(0.4)] 2 / Co(0.4)/Ta(0.2)/CoFeB(1)/ MgO(1.1)/CoFeB(1.8)/Ta(5)/ Ru (5), where the numbers are in nanometers (FIG.1a).Films are deposited at room temperature by dc/rf magnetron sputtering on a thermally oxidized Si substrate.The devices are fabricated into a circular shape with a 40-80 nm diameter using electron beam lithography and Ar ion milling and annealed at 300 • C for 1 hour by applying a 0.4 T magnetic field in the perpendicular direction.The average tunnel magnetoresistance ratio (TMR) and resistance area product (RA) are 65% and 4.7 Ωµm 2 , respectively.The discrete sMTJs used in this work are first cut out from the wafer and the electrode pads of the sMTJs are bonded with wires to IC sockets.The following parameters are measured by sweeping DC current to the sMTJ and measuring the voltage.The resistance of the P state R P is 4.4-5.7 kΩ, the resistance of the AP state R AP is 5.9-7.4kΩ, and the current at which P/AP fluctuations are 50% is defined as I 50/50 , in between 14-20 µA.At the output of the new p-bit design, we use an extra branch with a bipolar junction transistor (BJT) that acts as a buffer to the peripheral module (PMOD) pins of the Kintex UltraScale KU040 FPGA board.Given the electrostatic sensitivity of the sMTJs, this branch also protects the circuit from any transients that might originate from the FPGA.

Digital synthesis flow
HDL codes are converted to gate-level models using the Synopsys Design Compiler.Conversion from these models to Spice netlists is done using Calibre Verilog-to-LVS.Netlist post-processing is done by a custom Mathematica script to make it HSPICE compatible.Details of the synthesis flow (shown in FIG.4), followed by HSPICE simulation results for functional verification and power analysis are provided in Supplementary Sections XIII, XIV and XV.
where rand(-1, 1) represents a random number drawn from the uniform distribution in [−1,1], β is the inverse algorithmic temperature, and I i is the local field of p-bit "i" received from its neighbors.τ N in this equation is defined as the neuron evaluation time [65].For the typical choice of 2-local (Ising-like) energy functions, I i is given by: where J ij are the weights and h i is the bias term for each individual p-bit.τ S represents the synapse evaluation time.If the network of p-bits is symmetric (J ij = J ji ), it is possible to define the following energy function: In such a network of N p-bits, there are 2 N states, S = {1, 2, . . ., 2 N }, and each state j ∈ S is visited according to the Boltzmann-Gibbs distribution: The coupled evolution of Supplementary Eq. ( 1) and Supplementary Eq. ( 2) represents a dynamical system and a wide variety of problems can be mapped onto this system, including combinatorial optimization, machine learning and quantum simulation problems, all phrased in terms of powerful physics-inspired Monte Carlo algorithms [6,8].In this paper, we focus on the two settings of probabilistic inference and Boltzmann learning.In either case, we will be interested in a fixed β value (typically 1) and sample from the corresponding Boltzmann-Gibbs distribution of the network.We will show how the asynchronous and truly random bits generated by the sMTJs represent a low-level realization of physics-inspired probabilistic computation.

II. CHARACTERIZING SMTJS THROUGH P-BITS
We first characterize the statistics of the sMTJs (in Supplementary Figure 1) by making measurements on the p-bit circuit described later in Supplementary Figure 2b.The fluctuations of the p-bit circuit are controlled entirely by the sMTJs.We observe the outputs of 5 sMTJ-based p-bits and characterize their rate of fluctuations from event times and autocorrelations.In the main paper FIG.1g shows the fluctuations of p-bits #1-5 at I 50/50 , which is the current at which the sMTJs show 50/50 fluctuations for the high-and low-resistance states.The rate of fluctuations provides an estimate of the neuron evaluation time τ N of Supplementary Eq. ( 1) after which the p-bit produces a new and independent random bit.Even though the sMTJs were fabricated under the same conditions, slight differences in their volumes and shapes critically affect their relaxation times.The relaxation times τ P,AP obeys a Néel-Arrhenius law [40], described by τ = τ 0 exp(∆/k B T ), where τ 0 is the attempt time and ∆ is the energy barrier of the nanomagnet.In this paper, our magnets are not truly zerobarrier and their fluctuation rates exponentially depend on the energy barrier, ∆, however, detailed theoretical analysis shows that this dependence is much less pronounced in low barrier nanomagnets [22,23].Moreover, as we experimentally show in Sections 3-4, p-bits in symmetric networks are agnostic to update orders [3,65,66], therefore variations in these time scales should not play a prohibitive role in scaled implementations of p-computers (see Supplementary Section VIII for another experiment showing this update order invariance).
In Supplementary Figure 1b, we calculate the normalized autocorrelation of the p-bits using a 300 s time window with a sampling rate of 3.16 kHz, collecting N ≈ 949000 samples.The discrete autocorrelation function is defined as: where V [n] is the discrete sampled voltage read from the oscilloscope, and m represents the discretized lag time.Since autocorrelations typically show exponential decay [22], we extract an autocorrelation time τ ac , by measuring and fitting the autocorrelation to a continuous function C(t lag ) = exp(−t lag /τ ac ).As an additional measure to τ ac , we also determine the frequency of magnetic relaxations that are characterized by an event time.The event time, t event , is defined as the time between one event to the next.Measuring the distribution of event times from parallel and antiparallel configurations, we observe that the p-bit outputs are distributed according to the exponential distribution, indicating a Poisson process [25].Overall, p-bits show a good degree of variation in the relaxation and autocorrelation times as well as in their TMR values and their 50/50 currents.These values are summarized in Supplementary Table 1.

III. A VOLTAGE-COMPARATOR BASED P-BIT
In this section, we describe a new p-bit circuit we designed in this work, comparing its characteristics to an earlier design implemented in [28].Supplementary Figure 2a,b show both these designs.The earlier design was implemented in several small-scale experimental demonstrations using perpendicular sMTJs [28,29,67].In the original theoretical proposal [15], however, circular or elliptical in-plane nanomagnets were used.In-plane low barrier magnets are very hard to pin, requiring spin polarized currents of ≈100-500 µA or more for typical parameters [23].If the magnetization provides continuous randomness providing all resistance values between R P to R AP , this allows a faithful realization of Supplementary Eq. ( 1) as carefully discussed in [68].In such a case, spin-transfer-torque pinning is an unnecessary distraction, other than causing a read disturbance.Indeed, the fastest experimental p-bits are based on in-plane magnetic tunnel junctions [26,27] and variations of the design proposed in Ref. [15] may still be useful in future implementations.
The experimental demonstrations including the present work have so far been primarily focused on perpendicular sMTJs to keep fluctuation speeds slow for practical reasons.Perpendicular sMTJs are easily pinned with spin currents around ≈10 to 20 µA for typical parameters [68].Unlike in-plane sMTJs, however, perpendicular sMTJs switch telegraphically and they do not provide a uniform resistance between the two extremes.By a fortunate coincidence, the presence of the uncompensated dipolar field from the reference layer and easy pinning of perpendicular sMTJs appear to allow the realization of Supplementary Eq. (1) in hardware [28], since the spin-torque pinning changes the 50/50 fluctuations of the sMTJ [67].
In the design shown in Supplementary Figure 2a, the comparator has a fixed reference voltage.This means that as a function of V IN , the drain voltage swings between two extremes, V D,AP = V DD − R AP I MTJ and V D,P = V DD − R P I MTJ .As shown

B. RNG Unit
As shown in Algorithm 1 in the Supplementary information, we used three different types of RNGs: linear feedback shift register (LFSR), Xoshiro [47] and sMTJ-driven LFSR to compare the quality of randomness.We used 32-bit RNGs and compared the RNG outputs with the 32-bit LUTs to implement Supplementary Eq. ( 1).LFSR involves a linear shift operation on all the bits and an XNOR operation on some bits based on the selected taps.Unique seeds were used for each RNG and unique taps were used for RNGs in each p-bit block while ensuring maximal-length outputs.In contrast, a Xoshiro RNG involves linear shift and rotation operations on 32-bit words as well as XORing between subsequent words.In the all-digital versions, LFSR and Xoshiro RNGs are driven by digital clocks.However, in the hybrid design where sMTJs serve as the clocks for the RNGs, sMTJ + LFSR is considered as a new RNG unit.Using sMTJ-based p-bit removes the need for the 32-bit RNG and the 32-bit LUT.

C. Clocking and Sampling Unit
For the fully digital CMOS p-bits, each p-bit and its RNG is driven by system clocks generated on the low-voltage differential signaling (LVDS) clock-capable pins of the FPGA.These clocks are generated using Xilinx LogiCORE™ IP clocking wizard and mixed-mode clock manager (MMCM) module.Each of the five clocks operates at 15 MHz with shifted phases and is highly accurate with low jitter noise.In the Boltzmann learning example and for the NIST tests, to make these clocks comparable to sMTJs, we used frequency divider circuits to slow down the clocks to ≈ 2 kHz.For the sMTJ-driven p-bits, sMTJs replace the FPGA clocks and drive the p-bits and the RNGs externally using the peripheral module (PMOD) interface.FPGA samples the sMTJ outputs with a fast system clock of 75 MHz.

D. Data Programming and Acquisition Unit
We used MATLAB 2022a as the host program to read-write data to and from the FPGA through the USB-JTAG interface.MATLAB communicates with the FPGA board via AXI4 (Advanced eXtensible Interface 4) protocol where MATLAB works as the AXI master to drive a slave memory-mapped registered bank and Block RAMs (BRAM) inside the FPGA.We used airHDL [69], a memory management tool to assign the memory addresses for the register bank and the BRAMs.The weights J, h of the full adder circuit are programmed through MATLAB.In the inference and the Boltzmann learning example, the p-bit outputs were read from MATLAB as batches of sweeps that were initially sampled at 2kHz and stored in a BRAM.For the NIST tests, however, we sampled and stored all the sweeps in the BRAM at a much higher frequency (≈10 kHz) compared to the clock frequency of ≈ 2 kHz and then downsampled the data to a designated frequency.This procedure ensured that we did not lose any samples.MATLAB reads the BRAM data in burst mode and performs the downsampling.

VI. INFERENCE ON THE PROBABILISTIC FULL ADDER
Inside the FPGA, we construct digital p-bits that behave according to Supplementary Eq. ( 1) and interconnections between p-bits that behave according to Supplementary Eq. (2).Each p-bit has a PRNG, a LUT for the hyperbolic tangent function and 10-bit weights in fixed-precision, 1 sign, 6 integer, and 3 fractional bits (s [6][3]).Supplementary Eq. ( 2) is implemented by a multiply-accumulate unit inside the FPGA, whose multiplication reduces to simple multiplexing since a given weight J ij is either taken or not if m j is 0 or 1.An important consideration in ensuring the p-bit network reaches the equilibrium is the necessity of fast synapse times (τ s ) compared to neuron times τ n [65].In our context this requirement (τ s < τ n ) is naturally satisfied because the combinational logic inside the FPGA, which computes Supplementary Eq. ( 2) with about 10 ns delays is orders of magnitude faster than both our deliberately slowed digital clocks and our sMTJs.In scaled and integrated implementations with fast p-bits with GHz fluctuations, this necessity requires careful design.In the case sMTJ clocks, there is also the theoretical possibility of parallel updates by simultaneously switching sMTJs.Practically this is not a concern due to the extremely low probability of such an event which would be washed over thousands of samples anyway.Our results with sMTJs in FIG.2c,d and in FIG.3c,d indicate that the sMTJs reproduce the ideal distributions well.
A full block diagram of the FPGA unit is shown in FIG.3a.To rule out any spurious correlations, the starting states (seeds) used for the LFSR and Xoshiro are randomized for each p-bit.LFSRs also use unique sets of random taps while ensuring maximum-length outputs.
In this setting, for each RNG, we cumulatively sampled 10 6 states from the p-bits starting from a random initial state.A system state out of 5-p-bits can be defined from 0 to 2 N − 1 such that state 0 is p 1 p 2 p 3 p 4 p 5 = 00000 and state 31 is p 1 p 2 p 3 p 4 p 5 = 11111.We define the single update of each p-bit according to Supplementary Eq. ( 1)-(2) as a sweep.Due to their digital nature, defining exact times to perform a sweep for LFSR and Xoshiro is straightforward.With a driving clock frequency of 2 kHz, we perform one sweep and then record it as a new state.However, sampling states from sMTJ-clocked LFSR p-bits is not straightforward due to the analog nature of fluctuations of the sMTJs.The relaxation time of the slowest sMTJ is ≈20 ms (50Hz) that is 40× slower than the sampling frequency of 2 kHz.For this reason, we collect 40× more data points from sMTJ-based p-bits and downsampled them to obtain independent samples.The oversampling and subsequent downsampling of sMTJ data is due to our 40× faster readout process, not an inherent sMTJ limitation, and is implemented 5 sMTJ-based p-bits we characterized (Supplementary Figure 1) and sampled them at 2 kHz.To train the full adder on the LFSR and Xoshiro RNG-based p-bit network, we took 400 sweeps per epoch for a total of 500 epochs.For the sMTJ-clocked LFSR based p-bit network, we took 16000 sweeps per epoch and a total of 500 epochs.We took 40× more sweeps for the sMTJ because they were sampled at 2 kHz (40× faster than their autocorrelation) in order to produce the same number of independent sweeps for all solvers.
In Supplementary Figure 4 we provide the full histograms (32-states) of the full adder for the sampling and learning experiments.Only parts of the histograms are shown in the main text for clarity.In Supplementary Figure 8 we show that this bias is not random and never goes away, rather it is a consistent systematic bias by starting the LFSRs at different uniform random initial conditions (seeds) and different maximum-length taps.In Supplementary Figure 5, we perform additional experiments to study how undirected Boltzmann networks are invariant to update orders in the system, significantly easing experimental difficulties in scaled-out implementations of p-bit networks.In this experiment, a fully connected undirected p-bit network had nodes clocked asynchronously with independent sMTJs.Even with only 10 4 samples, a very close match is seen between the distributions corresponding to sMTJ-clocked LFSRs and Gibbs sampling in simulation.It is important to note that Gibbs sampling was performed with randperm updates, where a new update order (out of 4!=24 possibilities) was sampled (without replacement) at each iteration.Despite this highly random updating scheme, both Gibbs (ideal simulation) and the sMTJ-clocked system eventually reach the ground truth, represented by the Boltzmann distribution.This remarkable feature of update invariance in Boltzmann networks [66] significantly eases fabrication difficulties in eventual sMTJ-based large-scale p-bit networks.Going further, the investigation of "time-to-equilibrium" that determines the model averages and correlations in Algorithm 1 remains to be one of the future challenges of the field.

IX. SAMPLING FROM ARBITRARY DISTRIBUTIONS
In this section, we show how the CD algorithm can be used to create networks of p-bits that can sample from arbitrary distributions up to M-bit accuracy.We start by defining a PDF for a given distribution.As an example, we choose a normal distribution with µ = 16 and σ = 4 to be approximated by a N = 5 p-bit network with 2 N = 32 states.Once the PDF is specified, one approach is to create a truth table matrix of size V = N T × N , where the rows of the truth table are repeated according to the specified PDF.Assuming a fully visible network, we then calculate the data correlations from the truth table, V † V .We then perform the standard contrastive divergence algorithm.Supplementary Figure 6 shows inference and learning of a normal and exponential distribution with Mersenne Twister-based high-quality RNGs and LFSR-based low-quality RNGs.Even though the Mersenne-based learning looks marginally better, we do not see a strong LFSR bias in such low-dimensional

Bias
Supplementary Figure 8. Demonstration of consistent systematic bias in sampling the 5 p-bit full adder using multiple sets of seeds and taps for the LFSR.The seeds are chosen uniform randomly and the taps are chosen to provide the maximum length of the LFSR.Irrespective of the seeds and taps, the bias never goes away, indicating a systematic bias.
100 MHz to 1 GHz were tested, with 1 GHz being used for all results reported in this work.As seen from Supplementary Figure 10a,b, the experimental data obtained from HSPICE by varying the p-bit inputs (8-bit LUT input) of the synthesized circuit falls on the theoretically expected 1 + tanh(x)/2 for both LFSR and Xoshiro-based pbits.Varying the p-bit input allows us to tune the probability with which the p-bit fluctuates in accordance with the sigmoid seen in Supplementary Figure 10a,b.In Supplementary Figure 10c,d, we observe that across 25000 experimental samples, the distribution of outputs obtained from the digitally synthesized PRNG (LFSR or Xoshiro) are uniformly random.These two results verify the functionality of the digital p-bit synthesized using ASAP7 by a transistor-level simulation performed in HSPICE.

XIV. FUNCTIONAL VERIFICATION OF P-BITS AND PRNGS
fully-integrated and sMTJ-based computers with N = 10 6 p-bits can reach sampling throughputs of 10 6 flips/ns, 5 orders of magnitude faster than typical TPU/GPUs.TPU/GPU references discussed have been plotted on a power consumption versus sampling throughput scale in [77] and [74].
For sampling throughput projections at large scale, we use N/τ , assuming each p-bit flips independently of each other.This basic formula assumes that each flip is communicated to neighboring p-bits before a new flip is attempted, as otherwise flips may not be useful.If the network density scales as O(k • N ) with some small k corresponding to sparse networks, the fast communication assumption is warranted and ideal parallelism can be achieved.
For FPGA-based p-computers, the total power consumption is around 30W, including peripheral and unrelated circuits beyond the digital synthesis of our design.For heterogeneous computers of the type we consider in this work (5 MTJs + FPGA), the total power is also dominated by the FPGA power and is also around 30W. Unlike the present work where we used sMTJs to clock digital PRNGs such as LFSRs, in the future we envision the sMTJ-based analog p-bits as standalone blocks interacting through a CMOS underlayer without any seeding of CMOS PRNGs.In terms of power estimates for such fully sMTJ-based p-computers, detailed circuit simulations with experimentally established parameters indicate a power consumption of 10 µW per p-bit [68].For fully sMTJ-based p-computers with N = 10 6 , this would indicate a total p-bit power of 10 W, with an additional 10 W estimated synapse power [4,74].Considering how present day MRAM technology has been scaled up to 1 Gbit densities [78], integrating about 10 6 p-bits on top of CMOS should be reasonable.Supplementary Table 4. Benchmarking probabilistic hardware from device, circuit and system perspectives.]Benchmarkingprobabilistic hardware from device, circuit and system perspectives.For device comparisons, we focus on experimentally demonstrated sMTJ fluctuations.p-bits are circuits that use sMTJs to produce binary stochastic neurons with tunable probability with fluctuations at τ −1 rates.At the system level, we focus on the number of p-bits in a network (N) and sampling throughput, which is given by N/τ , for asynchronous systems with fast synapses computing Supplementary Eq. 2 (see text).Also at the system level, we report published data for GPUs/TPUs handling similar probabilistic sampling tasks.Projections are shown using ( †). (*) In heterogeneous computers of the type we consider in this paper, external sMTJ-based p-bits can drive a large number of digital p-bits in an FPGA or an ASIC.• RNG quality is deemed low / high for sampling problems rather than combinatorial optimization problems for which LFSR-based PRNGs seem sufficient [5].

Fig. 1 :
Fig. 1: Experimental setup for the CMOS + sMTJ probabilistic computer.(a) Stack structure of the stochastic magnetic tunnel junction (sMTJ).(b) The proposed sMTJ-based p-bit circuit with two branches whose outputs are provided to an operational amplifier.R ave is the average resistance of R P and R AP of the sMTJ.5 sMTJ-based p-bits provide tunable, truly random and asynchronous clocks to a digital field programmable gate array (FPGA).(c) Digital p-bits in the FPGA use lookup tables (LUT), comparators, synaptic weights, and pseudorandom number generators (PRNG).The clocks of the PRNG are driven by the truly random asynchronous outputs coming from the analog p-bits.(d) Pictorial representation of perpendicular sMTJ.(e) Image of a single p-bit circuit.(f) Image of the FPGA.The asynchronous clocks are input through the peripheral module (PMOD) pins.(g) Typical output of p-bits #1 to 5 using 5 sMTJs obtained from the p-bit circuit (see Supplementary Section III), showing variations in fluctuations.(h) Experimentally measured ⟨V OUT ⟩ the time average (over a period of 3 minutes) of the p-bit circuit output, as a function of DC input voltage V IN .The yellow squares are experimental data, and the blue dashed line is a fit of the form⟨V OUT ⟩ = 1/2 V ′ CC [tanh[β(V IN − V 0 )] + 1], where V 0 = 1.55 V, β = 3.43 V −1 , V ′ CC = 3 V is a reduced voltage from V CC = 5 V (see Supplementary Section IV).

Fig. 2 :
Fig. 2: Inference on a probabilistic full adder.(a) Fully-connected full adder network [39], where p-bits are clocked by the sMTJs.(b) Truth table of the full adder where Dec. represents the decimal representation of the state of [A B C in S C out ] from left to right.(c) Kullback-Leibler (KL) divergence between the ideal and measured distributions vs. the number of sweeps.Results are shown for LFSR-based p-bit (red line), sMTJ-clocked LFSR-based p-bit (blue line), and Xoshiro-based p-bit (green line).(d) Histogram for the measured and ideal distributions at the 10 6 sweep.The red, blue, and yellow bars show LFSR, sMTJ-clocked LFSR, and Boltzmann distribution, respectively.The histogram shows all 8 high probability states denoted in (b) and with a clear bias for the LFSR distribution (see Supplementary Section VII for full histograms for all PRNGs including Xoshiro).

FPGAFig. 3 :
Fig. 3: Learning deep Boltzmann machines.(a) The architecture of the p-computer for learning.The digital p-bits in FPGA are fed by sMTJ-based p-bits output similar to probabilistic inference.The weights J ij and biases h i are updated in the CPU for a specified number of epochs.(b) (Top) The 32-node Chimera graph is used as a deep BM. (Bottom) An asynchronous clocking scheme is shown with node coloring.(c) KL divergence as a function of the number of epochs for LFSR (red line), LFSR clocked by sMTJ-based p-bit (blue line), and Xoshiro (green line) (d) The distribution of full adder with learned weights and biases at epoch = 400 where the number of sweeps per epoch = 400 for LFSR-only and the number of sweeps per epoch = 16000 for sMTJ-clocked LFSR.The Boltzmann distribution was obtained with β = 3.The red, blue, and yellow bars show LFSR and LFSR clocked by sMTJ-based p-bit, and Boltzmann, respectively.The histogram shows 4 correct (0, 6, 18, 25) and 2 incorrect (2, 29) states, out of the 32 possible states.sMTJ-based p-bit closely approximates the ideal Boltzmann distribution whereas the LFSR underestimates correct states and completely fails with states 2 and 29 (see Supplementary Section VII for full histograms for all PRNGs including Xoshiro).
Supplementary Figure 1.sMTJ-based p-bit characterization.(a) Histogram of ln(MR) (number of magnetic relaxation (MR) events) as a function of the event times tevent for p-bit #2, following the exponential distribution 1/τP,AP exp(−tevent/τP,AP), expected from a Poisson process (red lines).(b) Autocorrelation of sMTJ-based p-bits #1-5.We define τac as the time at which the normalized autocorrelation decays to 1/e.

Supplementary Figure 4 .
Extended data for probabilistic sampling and learning experiments, showing all 32 states of the 5-bit full adder for LFSR and sMTJ-clocked LFSR.(a) Sampling data is taken at 10 6 sweeps, (b) Learning data is taken at epoch = 400 where the number of sweeps per epoch = 400 for LFSR-only and the number of sweeps per epoch = 16000 for the sMTJ-clocked LFSR.(c) Sampling data for Xoshiro (d) Learning data for Xoshiro VIII.UPDATE ORDER INVARIANCE OF BOLTZMANN NETWORKS

Supplementary Figure 10 .Supplementary Figure 11 .
(a) Verifying the functionality of 32-bit LFSR and (b) Xoshiro-based p-bits synthesized using ASAP7 PDK.Input current is converted from the 2's complement s43 representation to decimal equivalent.The y-axis indicates the probability of the p-bit being in a 1 state calculated over 2500 clock cycles of a transient p-bit simulation generating random bits.(c) Verifying the functionality of 32-bit LFSR and (d) Xoshiro PRNGs synthesized using ASAP7 PDK.32-bit binary outputs were mapped to their decimal equivalent, and a 25000 clock cycle transient simulation generated random 32-bit values where f clk = 1 GHz.Input profile to activate the LUT: 1 GHz positive and negative pulse are applied to the input of the p-bit to activate LUT transistors.Insets show representative transient simulations for instantaneous power consumption for LFSR and Xoshiro-based p-bits.

Table 1 .
Table of results for all 5 p-bits reporting mean relaxation time τ = √ τPτAP, τac, TMR and I 50/50 of sMTJs.I 50/50 is defined as the absolute value of the current at which the perpendicular sMTJs show 50/50 fluctuations by canceling the uncompensated dipolar field resulting from the fixed layer.