Probabilistic computing with p-bits

Jan Kaiser1, a) and Supriyo Datta1

Elmore Family School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN,
47906 USA
(Dated: 14 October 2021)
Digital computers store information in the form of bits that can take on one of two values 0 and 1, while quan-
tum computers are based on qubits that are described by a complex wavefunction whose squared magnitude
gives the probability of measuring either a 0 or a 1. Here we make the case for a probabilistic computer based
on p-bits which take on values 0 and 1 with controlled probabilities and can be implemented with specialized
compact energy-efficient hardware. We propose a generic architecture for such p-computers and emulate
systems with thousands of p-bits to show that they can significantly accelerate randomized algorithms used
in a wide variety of applications including but not limited to Bayesian networks, optimization, Ising models
and quantum Monte Carlo.
arXiv:2108.09836v2 [[Link]] 12 Oct 2021

I. INTRODUCTION Ideally the Kernel and data collector are pipelined so

that they can continually accept new random numbers
Feynman1 famously remarked “Nature isn’t classical, from the RN G4 , which is assumed to be fast and avail-
dammit, and if you want to make a simulation of nature, able in numerous numbers. The p-computer can then
you’d better make it quantum mechanical”. In the same provide Np fc samples per second, Np being the number
spirit we could say “Many real life problems are not deter- of parallel units9 , and fc the clock frequency. We argue
ministic, and if you want to simulate them, you’d better that even with Np = 1, this throughput is well in excess
make it probabilistic”. But there is a difference. Quan- of what is achieved with standard implementations on
tum algorithms require quantum hardware and this has either CPU (central processing unit) or GPU (graphics
motivated a worldwide effort to develop a new appropri- processing unit) for a broad range of applications and al-
ate technology. By contrast probabilistic algorithms can gorithms including but not limited to those targeted by
be and are implemented on existing deterministic hard- modern digital annealers or Ising solvers 8,10–17 . Interest-
ware using pseudo RNG’s (random number generators). ingly, a p-computer also provides a conceptual bridge to
Monte Carlo algorithms represent one of the top ten al- quantum computing, sharing many characteristics that
gorithms of the 20th century2 and are used in a broad we associate with the latter18 . Indeed it can implement
range of problems including Bayesian learning, protein algorithms intended for quantum computers, though the
folding, optimization, stock option pricing, cryptography effectiveness of quantum Monte Carlo depends strongly
just to name a few. So why do we need a p-computer ? on the extent of the so-called sign problem specific to the
A key element in a Monte Carlo algorithm is the RNG algorithm and our ability to ‘tame’ it19 .
which requires thousands of transistors to implement
with deterministic elements, thus encouraging the use of
architectures that time share a few RNG’s. Our work II. IMPLEMENTATION
has shown the possibility of high quality true RNG’s us-
ing just three transistors3 , prompting us to explore a dif- Of the three elements in Fig. 1, two are deterministic.
ferent architecture that makes use of large numbers of The Kernel is problem-specific ranging from simple op-
controlled-RNG’s or p-bits. Fig. 1 (a)4 shows a generic erations like addition or multiplication to more elaborate
vision for a probabilistic or a p-computer having two operations that could justify special purpose chiplets20 .
primary components: an N-bit random number gener- Matrix multiplication for example could be implemented
ator (RNG) that generates N-bit samples and a Kernel using analog options like resistive crossbars16,21–23 . The
that performs deterministic operations on them. Note data collector typically involves addition and could be
that each RNG-Kernel unit could include multiple RNG- implemented with adder trees. The third element is prob-
Kernel sub-units (not shown) for problems that can ben- abilistic, namely the N-bit RNG which is a collection of
efit from it. These sub-units could be connected in series N 1-bit RNG’s or p–bits. The behavior of each p–bit can
as in Bayesian networks (Fig. 2 (a) or in parallel as done be described by24
in parallel tempering5,6 or for problems that allow graph
coloring7 . The parallel RNG-Kernel units shown in Fig. si = Θ [ σ(Ii − r) ] (1)
1 (a) are intended to perform easily parallelizable oper-
ations like ensemble sums using a data collector unit to where si is the binary p-bit output, Θ is the step function,
combine all outputs into a single consolidated output. σ is the sigmoid function, Ii is the input to the p-bit
and r is a uniform random number between 0 and 1.
Eq. (1) is illustrated in Fig. 1 (b). While the p-bit
output is always binary, the p-bit input Ii influences the
a) kaiser32@[Link]
mean of the output sequence. With Ii = 0, the output

1
 (a) Architecture for a p-computer (b) 1-bit RNG: p-bit

N-bit
Kernel
RNG Characteristic
LFSR
+LUT
p-bit output
tanh(Input)

Data
Moving Average

Implementation
Collector

Output
N-bit 3T +
Kernel Input s-MTJ
RNG

FIG. 1. Probabilistic computer: (a) Overall architecture combining a probabilistic element (N-bit RNG) with deterministic
elements (kernel and data collector). The N-bit RNG block is a collection of N 1-bit RNG’s, or p-bits. (b) p-bit: Desired
input-output characteristic along with two possible implementations, one with CMOS technology using linear feedback shift
registers (LFSR’s) and lookup tables (LUT’s)8 and the other using three transistors and a stochastic magnetic tunnel junction
(s-MTJ)3 . The first is used to obtain all the results presented here, while the second is a nascent technology with many
unanswered questions.

is distributed 50 − 50 between 0 and 1 and this may be like s-MTJ’s, naturally generate true random numbers
adequate for many algorithms. But in general a non- with infinite repetition period and ideally require only 3
zero Ii determined by the current sample is necessary to transistors and 1 MTJ.
generate desired probability distributions from the N-bit A simple performance metric for p-computers is the
RNG-block. ideal sampling rate Np fc mentioned above. The re-
One promising implementation of a p-bit is based on a sults presented here were all obtained with an field-
stochastic magnetic tunnel junction (s-MTJ) as shown in programmable gate array (FPGA) running on a 125 MHz
Fig. 1 (b) whose resistance state fluctuates due to ther- clock, for which 1/fc = 8 ns, which could be significantly
mal noise. It is placed in series with a transistor, and the shorter (even ∼ 0.1 ns25 ) if implemented with s-MTJ’s.
drain voltage is thresholded by an inverter3 to obtain a Furthermore, s-MTJ’s are compact and energy-efficient,
random binary output bit whose average value can be allowing up to a factor of 100 larger Np for a given area
tuned through the gate voltage VIN . It has been shown and power budget. With an increase of fc and Np , a
both theoretically25,26 and experimentally27,28 that s- performance improvement by 2-3 orders of magnitude
MTJ -based p-bits can be designed to generate new ran- over the numbers presented here may be possible with
dom numbers in times ∼ nanoseconds. The same cir- s-MTJ’s or other physics-based hardware.
cuit could also be used with other fluctuating resistors29 , We should point out that such compact p-bit imple-
but one advantage of s-MTJ’s is that they can be built mentations are still in their infancy30 and many ques-
by modifying magnetoresistive random access memory tions remain. First is the inevitable variation in RNG
(MRAM) technology that has already reached gigabit characteristics that can be expected. Initial studies sug-
levels of integration30 . gest that it may be possible to train the Kernel to com-
Note, however, that the examples presented here all pensate for at least some of these variations32,33 . Second
use p-bits implemented with deterministic CMOS ele- is the quality of randomness, as measured by statisti-
ments or pseudo-RNG’s using linear feedback shift regis- cal quality tests which may require additional circuitry
ters (LFSR’s) combined with lookup tables (LUT’s) and as discussed for example in Ref.27 . Certain applications
thresholding elements8 as shown in Fig. 1 (b). Such ran- like simple integration (Section III A) may not need high
dom numbers are not truly random, but have a period quality random numbers, while others like Bayesian cor-
that is longer than the time range of interest. The longer relations (Section III B) or Metropolis-Hastings methods
the period, the more registers are needed to implement that require a proposal distribution (Section III C) may
it. Typically a p-bit requires ∼ 1000 transistors30 , the have more stringent requirements. Third is the possible
actual number depending on the quality of the pseudo difficulty associated with reading sub-nanosecond fluctu-
RNG that is desired. Thirty-two stage LFSR’s require ations in the output and communicating them faithfully.
∼ 1200 transistors, while a Xoshiro128+31 would require Finally, we note that the input to a p-bit is an analog
around four times as many. Physics-based approaches, quantity requiring Digital-to-Analog Converters (DAC’s)

2
 (a) Bayesian Network (b) Correlations
F2 G61
100

GF M2 G62

GM G63 10−1
F1

|Correlation|
Correlation
M1 G64
10−2
RNG Kernel RNG
LUT LUT

Multiply & Accumulate

Self 5. generation
> >
Child 6. generation
LFSR LFSR 10 −3
CPT LUT

Grandchild Stranger
√
3. generation ∝ 1/ NS
LUT LUT 4. generation
> >
LFSR LFSR 101 102 103 104 105 106
to data collector to data collector Number of Samples,
NS NS

FIG. 2. Bayesian network for genetic relatedness mapped to a p-computer (a) with each node represented by one p-bit.
With increasing NS , the correlations (b) between different nodes is obtained more accurately.

unless the kernel itself is implemented with analog com- with N nodes can be mapped to a N-bit RNG-block feed-
ponents. ing into a Kernel which stores the conditional probabil-
ity table (CPT) relating it to the next generation. The
correlation between different nodes in the network can be
III. APPLICATIONS directly measured and an average over the samples com-
puted to yield the correct genetic correlation as shown
A. Simple integration
in Fig. 2 (b). Nodes separated by p generations have a
correlation of 1/2p . The measured absolute √ correlation
between strangers goes down to zero as 1/ Ns .
A variety of problems such as high dimensional integra-
tion can be viewed as the evaluation of a sum over a very This is characteristic of Monte Carlo algorithms,
large number N of terms. The basic idea of the Monte namely, to obtain results with accuracy ε we need Ns =
Carlo method is to estimate the desired sum from a lim- 1/ε2 samples. The p-computer allows us to collect sam-
ited number Ns of samples drawn from configurations α ples at the rate of Np fc = 125 MSamples per second if
generated with probability qα : Np = 1 and fc = 125 MHz. This is about two orders
of magnitude faster than what we get running the same
N NS algorithm on a Intel Xeon CPU.
X 1 X mα
M= mα ≈ (2) How does it compare to deterministic algorithms run
NS α=1 qα
α=1 on CP U ? As Feynman noted in his seminal paper1 , de-
terministic algorithms for problems of this type are very
The distribution {q} can be uniform or could be clev- inefficient compared to probabilistic ones because of the
erly chosen to minimize the standard deviation of the need to integrate over all the unobserved nodes {xB } in
34
estimate
√ . In any case the standard deviation goes down order to calculate a property related to nodes {xA }
as 1/ Ns and all such applications could benefit from a
p-computer to accelerate the collection of samples. Z
PA (xA ) = dxB P (xA , xB ) (3)

B. Bayesian Network
By contrast, a p-computer can ignore all the irrelevant
A little more complicated application of a p-computer nodes {xb } and simply look at the relevant nodes {xA }.
is to problems where random numbers are generated not We used the example of genetic correlations because it
according to a fixed distribution, but by a distribution de- is easy to relate to. But it is representative of a wide
termined by the outputs from a previous set of RN G0 s. class of everyday problems involving nodes with one-
Consider for example the question of genetic relatedness way causal relationships extending from ‘parent’ nodes
in a family tree35,36 with each layer representing one gen- to ‘child’ nodes37–39 , all of which could benefit from a
eration. Each generation in the network in Fig. 2 (a) p-computer.

3
 (a) Knapsack problem (b) Time to solution
CPU (DP) GPU (DP)
103 CPU (MCMC) p-comp. emul. (MCMC)

RNG Kernel CPU (Pisinger) p-comp. proj. (MCMC)

(s) (s)
LUT 101

solution
> ΔWeight Weight to

tosolution
Calculation < Capacity data
LFSR collector
Accept/ 10−1
LUT Reject

Time to
> ΔValue Value

Time
LFSR Calculation Function 10−3

10−5

102 103 104

Problem / nn
size,
Problem size

FIG. 3. Example of MCMC - Knapsack problem: (a) Mapping to the general p-computer framework. (b) Performance
of the p-computer compared to CPU implementation of the same probabilistic algorithm and along with two well-known
deterministic approaches. A deterministic algorithm like Pisinger’s40 which is optimized specifically for the Knapsack problem
can outperform MCMC. But for a given MCMC algorithm, the p-computer provides orders of magnitude improvement over
the standard CPU implementation.

C. Knapsack Problem with CPU (Intel Xeon @ 2.3GHz) and GPU (Tesla T4
@ 1.59GHz) implementations, using the probabilistic al-
gorithm. Also shown are two efficient deterministic algo-
Let us now look at a problem which requires random
rithms, one based on dynamic programming (DP), and
numbers to be generated with a probability determined
one due to Pisinger et al.40,43 .
by the outcome from the last sample generated by the
Note that the probabilistic algorithm (MCMC) gives
same RNG. Every RNG then requires feedback from the
solutions that are within 1% of the correct solution, while
very Kernel that processes its output. This belongs to
the deterministic algorithms give the correct solution.
the broad class of problems that are labeled as Markov
For the Knapsack problem getting a solution that is 99%
Chain Monte Carlo (MCMC). For an excellent summary
accurate should be sufficient for most real world applica-
and evaluation of MCMC sampling techniques we refer
tions. The p-computer provides orders of magnitude im-
the reader to Ref.41 .
provement over CPU implementation of the same MCMC
The knapsack is a textbook optimization problem de- algorithm. It is outperformed by the algorithm developed
scribed in terms of a set of items, m = 1, · · N , the mth , by Pisinger et al.40,43 , which is specifically optimized for
each containing a value vm and weighing wm . The prob- the Knapsack problem. However, we note that the p-
lem is to figure out which items to take (sm = 1) and computer projection in Fig. 3 (b) is based on utilizing
whichP to leave behind (sm = 0) such that the total value better hardware like s-MTJ’s but there is also significant
V = m vm sP m is a maximum, while keeping the total room for improvement of the p-computer by optimizing
weight W = m wm sm below a capacity C. We could the Metropolis algorithm used here and/or by adding
straightforwardly map it to the p–computer architecture parallel tempering5,6 .
(Fig. 1), using the RN G to propose solutions {s} at
random, and the Kernel to evaluate V, W and decide to
accept or reject. But this approach would take us to- D. Ising model
ward the solution far too slowly. It is better to propose
solutions intelligently looking at the previous accepted
Another widely used model for optimization within
proposal, and making only a small change to it. For our
MCMC is based on the concept of Boltzmann machines
examples we proposed a change of only two items each
(BM) defined by an energy function E from which one
time.
can calculate the synaptic function Ii
This intelligent proposal, however, requires feedback
from the kernel which can take multiple clock cycles. One Ii = β(E(si = 0) − E(si = 1)), (4)
could wait between proposals, but the solution is faster if
instead we continue to make proposals every clock cycle that can be used to guide the sample generation from
in the spirit of what is referred to as multiple-try Metropo- each RNG ‘i’ in sequence44 according to Eq. (1).
lis 42 . The results are shown in Fig. 34 and compared Alternatively the sample generation from each RNG

4
 (a) Transverse Field Ising Model (b) Correlations

10 250

RNG Kernel
LUT to

Multiply & Accumulate

> data

Weight LUT
Sample average
LFSR collector

LUT
>
LFSR

FIG. 4. Example of Quantum Monte Carlo (QMC) - Transverse Field Ising model: (a) Mapping to the general
p-computer framework. (b) Solving the transverse Ising model for quantum annealing. Subfigure (b) is adapted from B.
Sutton, R. Faria, L. A. Ghantasala, R. Jaiswal, K. Y. Camsari and S. Datta, IEEE Access, vol. 8, pp. 157238-157252, 2020
licensed under a Creative Commons Attribution (CC BY) license8 .

can be fixed and the synaptic function used to de- implementation would provide invertible logic that not
cide whether to accept or reject it within a Metropolis- only provides the output for a given input, but also gen-
Hastings framework45 . Either way, samples will be gen- erates all possible inputs corresponding to a specified
erated with probabilities Pα ∼ exp(−βEα ). We can solve output24,46,47 .
optimization problems by identifying E with the negative
of the cost function that we are seeking to minimize. Us-
ing a large β we can ensure that the probability is nearly E. Quantum Monte Carlo
1 for the configuration with the minimum value of E.
In principle, the energy function is arbitrary, but much Finally let us briefly describe the feasibility of using
of the work is based on quadratic energy functions defined p-computers to emulate quantum or q-computers. A q-
by a connection matrix Wij and a bias vector hi (see for computer is based on qubits that are neither 0 or 1, but
example8,10–17 ): are described by a complex wavefunction whose squared
X X magnitude gives the probability of measuring either a 0
E=− Wij si sj − hi si , (5) or a 1. The state of an n–qubit computer is described by
ij i a wavefunction {ψ} with 2n complex components, one
for each possible configuration of the n qubits.
ForPthis quadratic
energy function, Eq. (4) gives Ii = In gate-based quantum computing (GQC) a set of
β j Wij sj + hi , so that the Kernel has to per- qubits is placed in a known state at time t, operated on
form a multiply and accumulate operation as shown in with d quantum gates to manipulate the wavefunction
Fig. 4 (a). We refer the reader to Sutton et al.8 for through unitary transformations [U (i) ]
an example of the max-cut optimization problem on a
two-dimensional 90 × 90 array implemented with a p- {ψ(t + d)} = [U (d) ] · · · ·[U (1) ]{ψ(t)} (GQC) (6)
computer.
Eq. (4), however, is more generally applicable even if and measurements are made to obtain results with proba-
the energy expression is more complicated, or given by bilities given by the squared magnitudes of the final wave-
a table. The Kernel can be modified accordingly. For functions. From the rules of matrix multiplication, the
an example of a energy function with fourth order terms final wavefunction can be written as a sum over a very
implemented on an eight bit p-computer, we refer the large number of terms:
reader to Borders et al.30 . X (d)
A wide variety of problems can be mapped onto the (1)
ψm (t + d) = Um,i · · · · Uj,k ψk (t) (7)
BM with an appropriate choice of the energy function. i,··j,k
For example, we could generate samples from a desired
probability distribution P , by choosing βE = −`nP . An- Conceptually we could represent a system of n qubits and
other example is the implementation of logic gates by d gates with a system of (n × d) p-bits with 2nd states
defining E to be zero for all {s} that belong to the truth which label the 2nd terms in the summation in Eq.(7)48 .
table, and have some positive value for those that do Each of these terms is often referred to as a Feynman
not24 . Unlike standard digital logic, such a BM-based path and what we want is the sum of the amplitudes of

5
all such paths: Finally we note that quantum Monte Carlo methods,
both GQC and AQC, involve selective summing of Feyn-
nd
2
X man paths to evaluate matrix products. As such we
ψm (t + d) = A(α)
m (8) might expect conceptual overlap with the very active field
α=1 of randomized algorithms for linear algebra51,52 , though
the two fields seem very distinct at this time.
The essential idea of quantum Monte Carlo is to estimate
this enormous sum from a few suitably chosen samples,
not unlike the simple Monte Carlo stated earlier in Eq. IV. CONCLUDING REMARKS
(2). What makes it more difficult, however, is the so-
called sign problem 19 which can be understood intuitively In summary, we have presented a generic architecture
(α)
as follows. If all the quantities Am are positive then it is for a p-computer based on p-bits which take on values
relatively easy to estimate the sum from a few samples. 0 and 1 with controlled probabilities, and can be imple-
But if some are positive while some are negative with mented with specialized compact energy-efficient hard-
lots of cancellations, then many more samples will be ware. We emulate systems with thousands of p-bits to
required. The same is true if the quantities are complex show that they can significantly accelerate the implemen-
quantities that cancel each other. tation of randomized algorithms that are widely used for
The matrices U that appear in GQC are unitary with many applications53 . A few prototypical examples are
complex elements which often leads to significant cancel- presented such as Bayesian networks, optimization, Ising
lation of Feynman paths, except in special cases when models and quantum Monte Carlo.
there may be complete constructive interference. In gen-
eral this could make it necessary to use large numbers
of samples for accurate estimation. A noiseless quantum ACKNOWLEDGMENTS
computer would not have this problem, since qubits in-
tuitively perform the entire sum exactly and yield sam- The authors are grateful to Behtash Behin-Aein for
ples according to the squared magnitude of the result- helpful discussions and advice. We also thank Kerem
ing wavefunction. However, real world quantum com- Camsari and Shuvro Chowdhury for their feedback on
puters have noise and p-computers could be competitive the manuscript. The contents are based on the work done
for many problems. over the last 5-10 years in our group, some of which has
Adiabatic quantum computing (AQC) operates on been cited here, and it is a pleasure to acknowledge all
very different physical principles but its mathematical who have contributed to our understanding. This work
description can also be viewed as summing the Feynman was supported in part by ASCENT, one of six centers
paths representing the multiplication of r matrices: in JUMP, a Semiconductor Research Corporation (SRC)
program sponsored by DARPA.
[e−βH/r ] · · · ·[e−βH/r ] (AQC) (9)

This is based on the Suzuki-Trotter method described DATA AVAILABILITY

in Camsari et al.49 , where the number of replicas, r, is
The data that support the findings of this study are
chosen large enough to ensure that if H = H1 + H2 ,
available from the corresponding author upon reasonable
one can approximately write e−βH/r ≈ e−βH1 /r e−βH2 /r .
request.
The matrices e−βH/r in AQC are Hermitian, and their
elements can be all positive50 . A special class of Hamil-
tonians H having this property is called stoquastic and AUTHOR DECLARATIONS
(α)
there is no sign problem since the amplitudes Am in the
Feynman sum in Eq. (8) all have the same sign. Conflict of interests
An example of such a stoquastic Hamiltonian is the One of the authors (SD) has a financial interest in Lud-
transverse field Ising model (TFIM) commonly used for wig Computing.
quantum annealing where a transverse field which is
1 R.
quantum in nature is introduced and slowly reduced to P. Feynman, “Simulating physics with computers,” Int. J.
zero to recover the original classical problem. Fig. 4 Theor. Phys 21 (1982).
2 B. A. Cipra, “The Best of the 20th Century: Editors Name Top
adapted from Sutton et al.8 , shows a n = 250 qubit prob- 10 Algorithms,” SIAM News 33, 1–2 (2000).
lem mapped to a 2-D lattice of 250 × 10 = 2500 p-bits 3 K. Y. Camsari, S. Salahuddin, and S. Datta, “Implementing p-
using r = 10 replicas to calculated average correlations bits with Embedded MTJ,” IEEE Electron Device Letters 38,
between the z-directed spins on lattice sites separated by 1767–1770 (2017).
4 J. Kaiser, R. Jaiswal, B. Behin-Aein, and S. Datta, “Benchmark-
L. Very accurate results are obtained using Ns = 105
ing a Probabilistic Coprocessor,” arXiv:2109.14801 [cond-mat]
samples. However,these samples were suitably spaced to (2021), arXiv:2109.14801 [cond-mat].
ensure their independence, which is an important concern 5 R. Swendsen and J.-S. Wang, “Replica Monte Carlo Simulation

in problems involving feedback. of Spin-Glasses,” Physical review letters 57, 2607–2609 (1986).

6
 6 D. J. Earl and M. W. Deem, “Parallel tempering: Theory, ap- 23 M. Demler, “MYTHIC MULTIPLIES IN A FLASH,” , 3 (2018).
plications, and new perspectives,” Physical Chemistry Chemical 24 K. Y. Camsari, R. Faria, B. M. Sutton, and S. Datta, “Stochas-
Physics 7, 3910–3916 (2005). tic p -Bits for Invertible Logic,” Physical Review X 7 (2017),
7 G. G. Ko, Y. Chai, R. A. Rutenbar, D. Brooks, and G.-Y. Wei, 10.1103/PhysRevX.7.031014.
“Accelerating Bayesian Inference on Structured Graphs Using 25 J. Kaiser, A. Rustagi, K. Y. Camsari, J. Z. Sun, S. Datta,

Parallel Gibbs Sampling,” in 2019 29th International Conference and P. Upadhyaya, “Subnanosecond Fluctuations in Low-Barrier
on Field Programmable Logic and Applications (FPL) (IEEE, Nanomagnets,” Physical Review Applied 12, 054056 (2019).
Barcelona, Spain, 2019) pp. 159–165. 26 S. Kanai, K. Hayakawa, H. Ohno, and S. Fukami, “Theory of
8 B. Sutton, R. Faria, L. A. Ghantasala, R. Jaiswal, K. Y. Cam- relaxation time of stochastic nanomagnets,” Physical Review B
sari, and S. Datta, “Autonomous Probabilistic Coprocessing 103, 094423 (2021).
with Petaflips per Second,” IEEE Access , 1–1 (2020). 27 C. Safranski, J. Kaiser, P. Trouilloud, P. Hashemi, G. Hu, and
9 Here, we assume that every RNG-Kernel unit gives one sample J. Z. Sun, “Demonstration of Nanosecond Operation in Stochas-
per clock cycle. There could be cases where multiple samples tic Magnetic Tunnel Junctions,” Nano Letters 21, 2040–2045
could be extracted from one unit per clock cycle. (2021).
10 H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial opti- 28 K. Hayakawa, S. Kanai, T. Funatsu, J. Igarashi, B. Jinnai, W. A.

mization by simulating adiabatic bifurcations in nonlinear Hamil- Borders, H. Ohno, and S. Fukami, “Nanosecond Random Tele-
tonian systems,” Science Advances 5, eaav2372 (2019). graph Noise in In-Plane Magnetic Tunnel Junctions,” Physical
11 M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, Review Letters 126, 117202 (2021).
H. Tamura, and H. G. Katzgraber, “Physics-Inspired Optimiza- 29 O. Hassan, S. Datta, and K. Y. Camsari, “Quantitative Evalu-

tion for Quadratic Unconstrained Problems Using a Digital An- ation of Hardware Binary Stochastic Neurons,” Physical Review
nealer,” Frontiers in Physics 7 (2019), 10.3389/fphy.2019.00048. Applied 15, 064046 (2021).
12 K. Yamamoto, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, 30 W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari,

H. Teramoto, A. Sakai, S. Takamaeda-Yamazaki, and M. Mo- H. Ohno, and S. Datta, “Integer factorization using stochastic
tomura, “7.3 STATICA: A 512-Spin 0.25M-Weight Full-Digital magnetic tunnel junctions,” Nature 573, 390–393 (2019).
Annealing Processor with a Near-Memory All-Spin-Updates-at- 31 S. Vigna, “Further scramblings of Marsaglia’s xorshift genera-

Once Architecture for Combinatorial Optimization with Com- tors,” Journal of Computational and Applied Mathematics 315,
plete Spin-Spin Interactions,” in 2020 IEEE International Solid- 175–181 (2017).
State Circuits Conference - (ISSCC) (IEEE, San Francisco, CA, 32 J. Kaiser, R. Faria, K. Y. Camsari, and S. Datta, “Proba-
USA, 2020) pp. 138–140. bilistic Circuits for Autonomous Learning: A simulation study,”
13 M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, Frontiers in Computational Neuroscience 14 (2020), 10.3389/fn-
and H. Mizuno, “24.3 20k-spin Ising chip for combinational op- com.2020.00014.
timization problem with CMOS annealing,” in 2015 IEEE In- 33 J. Kaiser, W. A. Borders, K. Y. Camsari, S. Fukami,

ternational Solid-State Circuits Conference - (ISSCC) Digest of H. Ohno, and S. Datta, “Hardware-aware in-situ Boltzmann
Technical Papers (2015) pp. 1–3. machine learning using stochastic magnetic tunnel junctions,”
14 I. Ahmed, P.-W. Chiu, and C. H. Kim, “A Probabilistic Self- arXiv:2102.05137 [cond-mat] (2021), arXiv:2102.05137 [cond-
Annealing Compute Fabric Based on 560 Hexagonally Coupled mat].
Ring Oscillators for Solving Combinatorial Optimization Prob- 34 C. M. Bishop, Pattern Recognition and Machine Learning: All

lems,” in 2020 IEEE Symposium on VLSI Circuits (IEEE, Hon- ”Just the Facts 101” Material (Springer (India) Private Limited,
olulu, HI, USA, 2020) pp. 1–2. 2013).
15 S. Patel, L. Chen, P. Canoza, and S. Salahuddin, “Ising Model 35 R. Faria, J. Kaiser, K. Y. Camsari, and S. Datta, “Hardware

Optimization Problems on a FPGA Accelerated Restricted Design for Autonomous Bayesian Networks,” Frontiers in Com-
Boltzmann Machine,” arXiv:2008.04436 [physics] (2020), putational Neuroscience 15 (2021), 10.3389/fncom.2021.584797.
arXiv:2008.04436 [physics]. 36 R. Faria, K. Y. Camsari, and S. Datta, “Implementing Bayesian
16 F. Cai, S. Kumar, T. Van Vaerenbergh, X. Sheng, R. Liu, C. Li, networks with embedded stochastic MRAM,” AIP Advances 8,
Z. Liu, M. Foltin, S. Yu, Q. Xia, J. J. Yang, R. Beausoleil, W. D. 045101 (2018).
Lu, and J. P. Strachan, “Power-efficient combinatorial optimiza- 37 D. Koller and N. Friedman, Probabilistic Graphical Models:

tion using intrinsic noise in memristor Hopfield neural networks,” Principles and Techniques (MIT Press, 2009).
Nature Electronics 3, 409–418 (2020). 38 B. Behin-Aein, V. Diep, and S. Datta, “A building block for
17 S. Dutta, A. Khanna, A. S. Assoa, H. Paik, D. G. Schlom, hardware belief networks,” Scientific Reports 6, 29893 (2016).
Z. Toroczkai, A. Raychowdhury, and S. Datta, “An Ising Hamil- 39 K. Yang, A. Malhotra, S. Lu, and A. Sengupta, “All-Spin
tonian solver based on coupled stochastic phase-transition nano- Bayesian Neural Networks,” IEEE Transactions on Electron De-
oscillators,” Nature Electronics 4, 502–512 (2021). vices 67, 1340–1347 (2020).
18 K. Camsari and S. Datta, “Dialogue Concerning the Two Chief 40 S. Martello, D. Pisinger, and P. Toth, “Dynamic Programming

Computing Systems: Imagine yourself on a flight talking to an and Strong Bounds for the 0-1 Knapsack Problem,” Management
engineer about a scheme that straddles classical and quantum,” Science 45, 414–424 (1999).
IEEE Spectrum 58, 30–35 (2021). 41 X. Zhang, R. Bashizade, Y. Wang, S. Mukherjee, and A. R.
19 M. Troyer and U.-J. Wiese, “Computational Complexity and Lebeck, “Statistical robustness of Markov chain Monte Carlo ac-
Fundamental Limitations to Fermionic Quantum Monte Carlo celerators,” in Proceedings of the 26th ACM International Con-
Simulations,” Physical Review Letters 94, 170201 (2005). ference on Architectural Support for Programming Languages
20 T. Li, J. Hou, J. Yan, R. Liu, H. Yang, and Z. Sun, “Chiplet and Operating Systems (ACM, Virtual USA, 2021) pp. 959–974.
Heterogeneous Integration Technology—Status and Challenges,” 42 J. S. Liu, F. Liang, and W. H. Wong, “The Multiple-Try Method

Electronics 9, 670 (2020). and Local Optimization in Metropolis Sampling,” Journal of the
21 C. Liu, B. Yan, C. Yang, L. Song, Z. Li, B. Liu, Y. Chen, H. Li, American Statistical Association 95, 121–134 (2000).
Q. Wu, and H. Jiang, “A spiking neuromorphic design with re- 43 H. Kellerer, U. Pferschy, and D. Pisinger, Knapsack Problems

sistive crossbar,” in 2015 52nd ACM/EDAC/IEEE Design Au- (2004).

tomation Conference (DAC) (IEEE, 2015) pp. 1–6. 44 S. Geman and D. Geman, “Stochastic Relaxation, Gibbs Distri-
22 B. Li, Y. Shan, M. Hu, Y. Wang, Y. Chen, and H. Yang, butions, and the Bayesian Restoration of Images,” IEEE Trans-
“Memristor-based approximated computation,” in International actions on Pattern Analysis and Machine Intelligence PAMI-6,
Symposium on Low Power Electronics and Design (ISLPED) 721–741 (1984).
(IEEE, 2013) pp. 242–247.

7
45 W. K. Hastings, “Monte Carlo sampling methods using Markov Information & Computation 8, 361–385 (2008).
chains and their applications,” Biometrika 57, 97–109 (1970). 51 P. Drineas, R. Kannan, and M. W. Mahoney, “Fast Monte
46 Y. Lv, R. P. Bloom, and J.-P. Wang, “Experimental Demonstra- Carlo Algorithms for Matrices I: Approximating Matrix Multi-
tion of Probabilistic Spin Logic by Magnetic Tunnel Junctions,” plication,” SIAM Journal on Computing 36, 132–157 (2006).
IEEE Magnetics Letters 10, 1–5 (2019). 52 P. Drineas, R. Kannan, and M. W. Mahoney, “Fast Monte Carlo
47 N. A. Aadit, A. Grimaldi, M. Carpentieri, L. Theogarajan, Algorithms for Matrices II: Computing a Low-Rank Approxima-
G. Finocchio, and K. Y. Camsari, “Computing with Invert- tion to a Matrix,” SIAM Journal on Computing 36, 158–183
ible Logic: Combinatorial Optimization with Probabilistic Bits,” (2006).
IEEE International Electron Devices Meeting (IEDM) (To ap- 53 A. Buluc, T. G. Kolda, S. M. Wild, M. Anitescu, A. DeGen-

pear) (2021). naro, J. Jakeman, C. Kamath, Ramakrishnan, Kannan, M. E.

48 S. Chowdhury, K. Y. Camsari, and S. Datta, “Em- Lopes, P.-G. Martinsson, K. Myers, J. Nelson, J. M. Restrepo,
ulating Quantum Interference with Generalized Ising Ma- C. Seshadhri, D. Vrabie, B. Wohlberg, S. J. Wright, C. Yang,
chines,” arXiv:2007.07379 [cond-mat, physics:quant-ph] (2020), and P. Zwart, “Randomized Algorithms for Scientific Comput-
arXiv:2007.07379 [cond-mat, physics:quant-ph]. ing (RASC),” arXiv:2104.11079 [cs] (2021), arXiv:2104.11079
49 K. Y. Camsari, S. Chowdhury, and S. Datta, “Scalable [cs].
Emulation of Sign-Problem–Free Hamiltonians with Room-
Temperature $p$-bits,” Phys. Rev. Applied 12, 034061 (2019).
50 S. Bravyi, D. P. Divincenzo, R. Oliveira, and B. M. Terhal, “The

complexity of stoquastic local Hamiltonian problems,” Quantum