# Probabilistic Neural Circuits leveraging AI-Enhanced Codesign for Random
Number Generation

Suma G. Cardwell1, Catherine D. Schuman3 J. Darby Smith1, Karan Patel3, Jaesuk
Kwon 2, Samuel Liu 2,
Christopher Allemang1, Shashank Misra1, Jean Anne Incorvia 2 and James B.
Aimone1
1Sandia National Laboratories, Albuquerque, NM, USA
2University of Texas, Austin
3University of Tennessee, Knoxville

Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Stochasticity is ubiquitous in the world around us. However, our predominant
computing paradigm is deterministic. Random number generation (RNG) can be a
computationally inefficient operation in this system especially for larger
workloads. Our work leverages the underlying physics of emerging devices to
develop probabilistic neural circuits for RNGs from a given distribution.
However, codesign for novel circuits and systems that leverage inherent device
stochasticity is a hard problem. This is mostly due to the large design space
and complexity of doing so. It requires concurrent input from multiple areas
in the design stack from algorithms, architectures, circuits, to devices. In
this paper, we present examples of optimal circuits developed leveraging AI-
enhanced codesign techniques using constraints from emerging devices and
algorithms. Our AI-enhanced codesign approach accelerated design and enabled
interactions between experts from different areas of the microelectronics
design stack including theory, algorithms, circuits, and devices. We
demonstrate optimal probabilistic neural circuits using magnetic tunnel
junction and tunnel diode devices that generate an RNG from a given
distribution.

###### Index Terms:

Probabilistic Devices, Probabilistic Neural Circuits, AI-Enhanced Codesign,
Codesign, Neuromorphic computing

Figure 1: Top: In scientific computing applications today, uniform random
numbers need to be converted to distributions of interest through expensive
rejection sampling or related techniques. Bottom: By using device
stochasticity and co-designed circuits, we can directly convert biased coin
flips to our distributions of interest, avoiding repeated sampling loops.

##  I Introduction

Neuromorphic computing is an emerging paradigm that promises to alleviate the
challenges faced by current classical computing approaches by emulating key
computational principles from the brain. Recent work has shown that neural
networks can be used to sample probabilistic graphs [1, 2]. However, today’s
neuromorphic platforms have at best modest pseudo-random number generator
(PRNG) capabilities. In contrast, every synapse in the brain exhibits
stochastic vesicle release. Neuroscientists have observed that stochasticity
at the synapse and circuit scales allows for both synaptic development and
circuit functional dynamics, phenomena that are crucial for higher-level
cognitive functions. Recent theoretical work has also demonstrated that neural
algorithms can leverage RNGs in parallel to provide added capabilities to
probabilistic algorithms while leveraging the energy and time advantages of
neuromorphic parallelism [3]. We seek to leverage stochasticity in computing
by exploiting the underlying physics of emerging random number generator (RNG)
devices to build probabilistic neural architectures. We will leverage
stochasticity in computing, by making stochasticity ubiquitous and, crucially,
making it useful.

In many scientific computing applications today, uniform pseudo-random numbers
are converted to a desired distribution. When examining a particular collider
physics simulation [4], on the order of 270K uniform pseudo-random numbers are
used to simulate a single event, and up to billions of events must be
simulated to understand experimental data [5]. The CPU time for PRNGs in such
simulations is on the order of $40-50\%$ of the total compute time [5].
Contributing and adding to this cost is the necessity to convert uniform
samples to distributions of interest: as PRNGs typically produce uniformly
distributed numbers, it is necessary to convert these samples into
distributions of interest using techniques such as rejection sampling or
Markov Chain Monte Carlo algorithms (Fig. 1, top).

Direct random number generation leveraging stochastic devices can promise
significant energy savings for such applications. In this work, we explore
drawing samples directly from a simple, non-uniform distribution using true
random number generating (TRNG) device models (Fig. 1, bottom). We choose a
simple distribution among four outcomes that cannot be represented by the
outcome of two independent coin tosses alone. Rather, we employ a
probabilistic mixing through a “Hidden Dependence” process that can be
represented as a probabilistic circuit. The hidden process randomly selects
between various coin distributions. Though our choice of distribution and
hidden process circuit may seem arbitrary, they are important as they
demonstrate a methodology to correlate output among independent tosses.
Leveraging correlation of output is critical in direct drawing from non-
trivial distributions. This work showcases our results from developing neural
probabilistic circuits using AI-enhanced codesign leveraging emerging
neuromorphic devices for random number generation from a given distribution as
illustrated in Fig. 1.

##  II AI-enhanced Codesign

It is not immediately clear how to design circuits and tune device parameters
with emerging devices. To enable the design of circuit topology and circuit
parameters, we use evolutionary optimization. There are several key reasons
that we leverage evolutionary optimization for this task. First, it can be
used to optimize multiple objectives; in this case, we may be interested in
how accurate the circuit performs the task, while simultaneously minimizing
latency and/or energy usage. Second, it can be used to evolve both the
topology of the circuit and the parameters simultaneously. In this work, we
leverage the Library of Evolutionary Algorithms in Python (LEAP) [6] to evolve
the parameters of a probabilistic circuit, though in future work, for both
topology and parameter optimization, we will leverage Evolutionary
Optimization for Neuromorphic Systems (EONS) [7, 8] for evolving both non-
neuromorphic and neuromorphic circuits of these devices. Third, evolutionary
algorithms can potentially be used to creatively discover new ways to design
circuits or leverage the underlying device characteristics. Evolutionary
algorithms have been used for decades to design analog circuits [9, 10]. They
have been found to have surprising creativity in the design of novel solutions
to a variety of problems [11]. In this work, we use LEAP’s traditional
evolutionary algorithm with a real representation, where the parameters to be
evolved represent parameters associated with the devices themselves (described
in detail in Section III).

##  III Devices

We evaluated two devices for random number generation, namely Tunnel Diodes
(TD) and Magnetic Tunnel Junctions (MTJs). We included functional models with
abstracted device characteristics to evaluate the probabilistic circuit
developed by the AI-enhanced codesign framework.

###  III-A Tunnel Diodes

Figure 2: Tunnel Diode characteristics. (a) Schematic diagram of a tunnel
diode, and transport curve taken from a discrete component tunnel diode
(MP1103 from M-Pulse Microwave). The regions of tunneling current and
thermionic emission are shown. (b) Pulsing the current near the peak (black
arrow in (a)) produces the stated probability of being found in either the low
voltage or high voltage branch. Error bars are associated with the number of
samples.

The venerable tunnel diode (TD) has historically been used in high-speed
analog applications, and is a great candidate for a practical nanoscale random
number generator. They can be integrated into a standard CMOS (Complementary
Metal-Oxide-Semiconductor) fabrication process, can be shrunk to nanoscale
dimensions [12], and have CMOS-compatible current densities and voltages [13].
As shown in Figure 2 (a), a TD consists of a strongly n-doped and p-doped
junction, and conducts either by tunneling through or by thermionic emission
over the narrow depletion region. Under current bias, the random occupancy of
charge traps in the depletion region determine whether the device conducts
through tunneling ($\sim$ 0.05 V) or thermionic emission ($\sim$ 0.45 V).
Interpreting tunneling as tails, and thermionic emission as heads, tuning the
TD can be accomplished by tweaking the bias current near the region of the
peak, where lower currents bias the TD towards tunneling, and higher currents
towards thermionic emission, as shown in Figure 2 (b). The resulting coin
flips have been shown to pass established tests of statistical quality [14].
While the data shown in Figure 2 was taken with a large discrete component
device at the MHz speed of our data acquisition, both scaled and high-speed
devices have been demonstrated in the literature. The device in Ref. [12] has
a dimension under 100 nm, and consumes a power of 13 nW for tunneling, and 38
nW for thermionic emission. The time to generate a flip can be determined by
the time taken to charge the junction capacitance (0.25 fF) twice - once to
reset it, and once to charge it for the coin flip - 1.2 ns. Assuming no
additional overhead from parasitics, this produces an energy consumption of
0.05 fW for a heads, and 0.02 fW for a tails. The device in Ref. [13] worked
at time scales as short as 50 ps, but is larger at 600 nm a side. The
correspondingly larger 0.8 mA current makes it simpler to integrate, with the
short operation time counteracting the larger power dissipation, limiting the
energy consumption to 50 fJ for heads and 20 fJ for tails. We use this last TD
device for resource estimation later in this manuscript.

###  III-B Magnetic Tunnel Junction

Figure 3: a) Schematic illustration of the VCMA-MTJ stochastic device
structure with various knobs to control the switching probability. T2-T4
applies a biased voltage on the MTJ structure. Spin transfer torque (STT)
(T1-T2) or spin Hall effect (SHE) (T1-T3) can bias the switching probability.
b) A and AP states are stable without applying a voltage. A biased voltage
eliminates stable states with a precessional switching to the middle in-plane
state. c) After removing the bias voltage, the device recovers to the stable
states with an unpredicted switching towards A or AP state.

Magnetic tunnel junctions (MTJs) have recently become a crucial component of
memory, in-memory computing, and probability bit (p-bit) device applications
[15, 16, 17]. An MTJ consists of an insulating tunnel barrier between two thin
ferromagnetic layers; due to spin-dependent electron transport, the MTJ has a
high and low resistance state depending on the P (parallel) and AP (anti-
parallel) orientation of the magnetization of the two ferromagnetic layers.
One of these layers is held at a fixed magnetization, while the other, the
free layer, can be switched via voltage, current, and heat. The thermally-
driven, stochastic nature of the switching of the free layer can lead to
generation of streams of random bits [18]. The MTJ can be set up as a
stochastic read device, by having a low-anisotropy magnetic free layer, near
the superparamagnetic limit, that randomly switches its magnetization at room
temperature. The random fluctuations of the MTJ resistance can then be read.
Alternatively, the MTJ can be set up as a stochastic write device, with a
stable-anisotropy free magnetic layer that has a probability of switching its
magnetization, and therefore the MTJ resistance, depending on the amplitude
and duration of an applied current pulse. In both of these configurations, the
weight of the probabilistic bit, i.e. how much time it spends in P or AP
states, can be controlled using constant applied fields and/or DC bias
voltages or currents. As an alternative to spin transfer torque (STT), the
spin Hall effect (SHE) can be used to change the state of the MTJ with low-
energy switching behavior as shown in Fig. 3 (a). We will refer to these
devices as MTJ-SHE.

An alternative way to generate stochastic bit streams with MTJs is using
voltage-controlled magnetic anisotropy (VCMA), where the anisotropy of the
free magnetic switching layer of the MTJ is modulated using voltage; this
effect has been extensively explored for memory applications. The VCMA effect
can exert unpredictable magnetization dynamics on the free layer at room
temperature, but it has not yet been fully explored for generating
controllable random bit streams. Here, we build numerical models of both the
VCMA-MTJs and the MTJ-SHE based on the Landau-Lifshitz-Gilbert (LLG) equation
[19, 20], modeling a standard perpendicular MTJ stack comprised of CoFeB (free
layer)/MgO/CoFeB (fixed layer), as shown in Fig. 3a. Fig. 3b-c depicts the
operation of the device using VCMA to generate a random bit stream: first, a
voltage is applied that reduces the perpendicular magnetic anisotropy and
sends the magnetization from out of plane to in-plane through precessional
oscillations (Fig. 3b). Then, the voltage is turned off, and when the
anisotropy pops back, the magnetization must choose one of two out-of-plane
directions to stabilize the MTJ in either a P or AP state (Fig. 3c). We will
refer to this device as MTJ-VCMA. Similarly, the MTJ-SHE uses the SHE to send
the magnetization to a middle-state when the SHE is on, and subsequently the
magnetization must recover to either a P or AP state.

The developed abstract models for both the MTJ-SHE and MTJ-VCMA capture the
devices functionality and energy consumption. These device models are utilized
when evaluating the fitness function of our evolutionary method and generated
bit stream from this analytical code is provided to the probabilistic circuit.

##  IV Probabilistic Circuits

An important application for probabilistic circuits is to generate RNGs from a
given distribution. There are many potential methods for taking a series of
two-state random outcomes, such as those produced by an MTJ or TD, and
converting them into a sample from a desired distribution. One approach would
be to utilize classic flipping methods [21]. Though, taking such approaches
may lead to tossing an unbounded number of arbitrarily tuned coins [21].
Another approach could come from identifying outcomes of a series of tossed to
coins to either a discrete distribution or a discretized continuous
distribution. Such an approach would guarantee a finite number of tosses
provided the distributions have finite support or are truncated in some
fashion. We will focus on sampling from discrete, finitely supported
distributions using such an identification. We will not review all such
methods of sampling nor argue that this is a superior one in any metric.
However, this test case demonstrates an important feature: we are able to
produce correlation among independent coin tosses to sample from a
distribution that cannot be sampled from two independent coin flips alone.

To mimic coin tossing, we consider the Bernoulli distribution. This
distribution is defined on two events, $0$(T) and $1$(H). We will say
$\mathbb{P}(1)=p$ and $\mathbb{P}(0)=1-p$. If we wish to sample from a
distribution with just two events, or a discretized continuous distribution
with two effective events, we need only make an identification of $H$ with
event 1 and $T$ with event 2 and assign the probability $p$ accordingly. Note,
for a discretized continuous distribution, this would mean assigning $p$ to be
the integral of the probability density function (pdf) over the bin
corresponding to event 1. If we have more than two events, more care may need
to be taken. For example, if we have four events, we may toss two coins in a
row. Let $p$ be the probability of heads for coin 1 and $q$ be the probability
of heads for coin 2. Since there are four distinct outcomes (HH, HT, TH, and
TT), we need only solve the appropriate system of equations to solve for the
probabilities $p$ and $q$. If you have three events, however, you will need to
arbitrarily select two of the four outcomes to represent a single one. In such
instances there may be many non-unique representations. This methodology
extends easily to any finite number of states.

However, situations may arise where there are no solutions for the coin heads
probabilities. In these situations, we will need to enforce correlation among
the outcomes of the coin tosses. There are many ways to correlate coin flips.
One way to correlate the output of coin flips among Bernoulli random variables
is through a hidden dependence model.

The “Hidden Dependence” model creates dependent random variables that are
connected by some underlying, possibly hidden, random variable as illustrated
in Fig. 4. This hidden variable can be thought of in many ways. It could, for
instance, represent some physical model. The dependence may induce correlation
among outputs. Relating this to COINFLIPS, the hidden variable could be the
physics of a device such as an MTJ or a TD. In our test case, the hidden
variable could be another coin whose outcome determines which set of two
additional coins are flipped. Hidden dependence coins can be flipped
simultaneously or in any order. This means from an outside viewpoint, the
coins appear to be independent tosses even though they represent dependent
random variables.

### Hidden Dependence Bernoulli Coins

In the “Hidden Dependence” scenario, there is a hidden process, stochastic or
deterministic, that controls the probability of heads among a collection of
coins, as shown in Fig. 4.

Figure 4: In the “Hidden Dependence” model, there is a hidden process,
stochastic or deterministic, that controls the probability of heads among a
collection of coins. In this cartoon, the hidden process chooses which set of
coins is flipped. The observer only sees a single set, the effective flipping
set.

The hidden process that controls the dependence among the effective visible
Bernoulli coins could be one of many options. It could be a deterministic
process that shuffles the used coin set; it could be a knob that tunes the
probability of heads among all coins; or it could be some more complicated
random determination. However, the key feature for this scenario is that the
effective visible set of coins can be flipped in any desired order (or all at
once). Any causation of output correlation and dependence is absorbed by
something a viewer does not see.

As an analogy, one could imagine a person that has two sets of coins in their
pocket, each set has a coin labeled coin 1 and a coin labeled coin 2. The
person encounters a series of individuals. When encountering an individual,
the person selects a set of coins and the individual is obliged to flip them.
If the person is strategic in how they choose the set of coins for each
individual to flip, they can induce a correlation in the effective output of
coin 1 and coin 2.

(a) $\omega_{1}$

(b) $\omega_{2}$

(c) $\omega_{3}$

Figure 5: Impact of multi-objective weights $\omega_{1}$, $\omega_{2}$,
$\omega_{3}$ on the KL divergence and energy usage of MTJ-SHE devices.

We now define a concrete example that we will develop a circuit for. Suppose
the distribution we would like to sample from is that of a four-sided die that
rolls $0$ with probability $1/2$, and rolls 1, 2, 3 with probability 1/6 each.
It’s easy enough to identify these outcomes with the outcomes of two coin
flips: 0 corresponds to HH, 1 to HT, 2 to TH, and 3 to TT. If $p$ is the
probability of $H$ for coin 1 and $q$ is the probability of H for coin 2, then
we need the following to hold:

$\displaystyle\begin{split}pq&=\displaystyle\frac{1}{2},\\\
p(1-q)&=\displaystyle\frac{1}{6},\\\ (1-p)q&=\displaystyle\frac{1}{6},\\\
(1-p)(1-q)&=\displaystyle\frac{1}{6}.\end{split}$ (1)

Despite the simple nature of this distribution, it is a quintessential example
of a system that requires the use of a hidden process. As hinted above, there
is a problem with this setup. This system of equations has no solution for $p$
or $q$. In order to realize the desired die roll, we will need to employ a
mechanism that alters the probabilities $p$ and $q$.

The set of equations (1) are ill defined as $p$ and $q$ are over-constrained.
To rectify this, we can relax the constraints by adding in additional
variables, namely by adding in a second set of coins. Returning to our
heuristic, we want our person handing out the sets of coins to hand out set 1
for a fraction $w$ of the time and set 2 for the remainder of the time. Set 1
has coin 1 with probability of heads $p_{1}$ and has coin 2 with probability
of heads $q_{1}$. Set 2 has probabilities of heads $p_{2}$ and $q_{2}$
respectively. This yields the following system of equations.

$\displaystyle\begin{split}wp_{1}q_{1}+(1-w)p_{2}q_{2}&=\displaystyle\frac{1}{2},\\\
wp_{1}(1-q_{1})+(1-w)p_{2}(1-q_{2})&=\displaystyle\frac{1}{6},\\\
w(1-p_{1})q_{1}+(1-w)(1-p_{2})q_{2}&=\displaystyle\frac{1}{6},\\\
w(1-p_{1})(1-q_{1})+(1-w)(1-p_{2})(1-q_{2})&=\displaystyle\frac{1}{6}\end{split}$
(2)

By altering between coin sets with probability $w$, we have induced a new
hidden process that controls the observed outcomes among the effective
observed coins. Notably, this introduction allows the system to be solved for
some choices of $w$. However, given that we have a choice in $w$, this means
we may not have a unique global solution for $p_{1}$, $p_{2}$, $q_{1}$, and
$q_{2}$.

While we concede that this exemplar is not tied to any specific application,
any practical distribution will not likely have outcome probabilities that may
be represented in a system similar to (1). This exemplar is meant to showcase
how one might correlate output and generate circuits to deal with such
problems. If the future of probabilistic computer design can use distributions
formed in such a manner, we will need the algorithmic expertise to handle
these types of over-constrained situations.

We model equation set (2) in Section V, and include the functional models of
the various emerging devices we studied in this work.

##  V AI-enhanced Probabilistic Circuits

We modeled the equations described in Section IV and leveraged LEAP to
generate the hidden weight ($w$) and probability weights ($p_{1}$, $q_{1}$,
$p_{2}$, and $q_{2}$) for the coinflip devices, where
$w,p_{1},q_{1},p_{2},q_{2}\in[0,1]$. For our examples we used LEAP’s
traditional evolutionary algorithm with a real representation, where the
parameters to be evolved represent parameters associated with emerging
devices, namely TD, MTJ-SHE, and MTJ-VCMA devices. We used tournament
selection, uniform crossover, and Gaussian mutation with a standard deviation
of $0.001$ and one expected mutation per individual. For our use-case this is
to find the optimal value of $w$, $p_{1}$,$q_{1}$, $p_{2}$, and $q_{2}$ such
that our desired objective is achieved. In this case, we have three key
objectives, described in our fitness function below.

Fitness function: The fitness function defines the objective of the
evolutionary algorithm. In this work, we leverage a multi-objective fitness
function, where we seek to have LEAP find the appropriate values for these
that result in minimizing the KL divergence value from the desired
distribution while simultaneously minimizing the difference between the
weights of coins and 0.5 (a fair coin), as well as minimizing energy usage
(based on the particular device). The KL divergence objective is specifically
measuring the difference from the given distribution. In particular, we
calculate the following:

$\displaystyle\begin{split}v(0)&=wp_{1}q_{1}+(1-w)p_{2}q_{2}\\\
v(1)&=wp_{1}(1-q_{1})+(1-w)p_{2}(1-q_{2})\\\
v(2)&=w(1-p_{1})q_{1}+(1-w)(1-p_{2})q_{2}\\\
v(3)&=w(1-p_{1})(1-q_{1})+(1-w)(1-p_{2})(1-q_{2})\\\ \end{split}$ (3)

We then set $p(0)=\frac{1}{2}$, $p(1)=\frac{1}{6}$, $p(2)=\frac{1}{6}$,
$p(3)=\frac{1}{6}$.

The KL divergence is then calculated:

$KL=\sum_{i=0}^{3}v(i)\log{\frac{v(i)}{p(i)}}$ (4)

For the given device model (MTJ-SHE, MTJ-VCMA or TD), we then calculate the
energy usage for devices with probabilities $p_{1}$, $p_{2}$, $q_{1}$, and
$q_{2}$. We sum these energy values and produce a single energy usage, $EN$.
To allow us to investigate the tradeoff between different objectives, we
include three objective weights $\omega_{1},\omega_{2},\omega_{3}$. Thus, our
overall fitness function is:

$\begin{split}f(w,p_{1},p_{2},q_{1},q_{2})&=\omega_{1}KL(w,p_{1},p_{2},q_{1},q_{2})\\\
&+\omega_{2}\left(\sum_{i=1}^{2}|p_{i}-0.5|+\sum_{i=1}^{2}|q_{i}-0.5|\right)\\\
&+\omega_{3}EN(p_{1},p_{2},q_{1},q_{2})\end{split}$ (5)

###  V-A Tradeoffs and Performance Evaluation

In Fig. 5, we can see the impact on the KL-divergence and the energy use for
different values for the three objective weights $\omega_{1}$, $\omega_{2}$,
$\omega_{3}$ for the MTJ-SHE device. Here, we can clearly see that by setting
the values for each weight differently, we can achieve different varying
levels of KL-divergence and energy usage. For example, a larger $\omega_{1}$
for the minimizing KL-divergence objective and a smaller $\omega_{3}$ for the
minimizing energy usage objective result in sets of device parameters that
achieve smaller KL-divergence values at the expense of using more energy. For
future co-design efforts, we can leverage this multi-objective approach to
tune the performance of our devices to achieve our desired tradeoffs in
performance and energy usage. We omit the plots for MTJ-VCMA and TD, where
there is less impact of the weight of the coins on the energy usage. In the
remainder of the paper, we focus on primarily minimizing KL divergence,
resulting in a focus on the following weight sets: $\omega_{1}=7500$,
$\omega_{2}=0.005$, and $\omega_{3}=0.5$.

###  V-B MTJ-VCMA Results

We selected the device configuration for MTJ-VCMA that gave the lowest KL-
divergence value. This device was optimized using $\omega_{1}=7500$,
$\omega_{2}=0.005$, and $\omega_{3}=0.5$. The resulting parameters found
through optimization over 1000 generations in LEAP for this device were
$w=0.233$, $p_{1}=0.237$, $p_{2}=0.781$, $q_{1}=0.199$, and $q_{2}=0.805$.
Fig. 6 plots the KL divergence and energy usage for the $n$ samples for this
MTJ-VCMA device across 10 different sets of samples. As seen Fig. 6,
increasing the number of samples lowers the KL-divergence from the desired
probability distribution. However, more samples come at the cost of increased
energy consumption. One of the resulting distributions with 2000 samples is
shown in Fig. 7.

Figure 6: KL divergence and energy usage vs. number of samples for the given
distribution with the MTJ-VCMA device. There were 10 trials conducted per
sample. Sample sizes investigated include 10, 50, 100, 200, 500, 1000, 1500
and 2000. Figure 7: Empirical distribution using MTJ-VCMA for 2000 samples in
a single run. The optimized values generated were: $w=0.306$, $p_{1}=0.448$,
$p_{2}=0.746$, $q_{1}=0.439$, and $q_{2}=0.749$.

###  V-C MTJ-SHE Results

We also selected the device configuration for MTJ-SHE that gave us the lowest
KL-divergence value. This device was optimized using $\omega_{1}=7500$,
$\omega_{2}=0.005$, and $\omega_{3}=0.5$. The resulting parameters found
through optimization over 1000 generations in LEAP for this device were
$w=0.306$, $p_{1}=0.448$, $p_{2}=0.746$, $q_{1}=0.439$, and $q_{2}=0.749$.
Fig. 8 plots the KL divergence and energy usage for the $n$ samples for this
MTJ-SHE device across 10 different sets of samples. As seen Fig. 8, increasing
the number of samples lowers the KL-divergence from the desired probability
distribution. However, more samples come at the cost of increased energy
consumption. One of the resulting distributions with 2000 samples is shown in
Fig. 9.

It is observed that there is an orders of magnitude energy reduction when SHE
is used compared to VCMA. In practice, SHE can require a DC external magnetic
field. In the model used, the SHE current is polarized perfectly in-plane;
experimentally, this is not always the case. The VCMA effect could provide
robustness at the cost of energy, since when VCMA is used, the magnetization
is dependent on anisotropy instead of spin current.

Figure 8: KL divergence and energy usage vs. number of samples for the given
distribution with the MTJ-SHE device. There were 10 trials conducted per
sample. Sample sizes investigated include 10, 50, 100, 200, 500, 1000, 1500
and 2000. Figure 9: Empirical distribution using MTJ-SHE for 2000 samples in a
single run. The optimized values generated were: $w=0.306$, $p_{1}=0.448$,
$p_{2}=0.746$, $q_{1}=0.439$, and $q_{2}=0.749$.

###  V-D TD Results

Again, we selected the device configuration for TD that gave us the lowest KL-
divergence value. This device was optimized using $\omega_{1}=7500$,
$\omega_{2}=0.005$, and $\omega_{3}=0.5$. The resulting parameters found
through optimization over 1000 generations in LEAP for this device were
$w=0.714$, $p_{1}=0.891$, $p_{2}=0.107$, $q_{1}=0.766$, and $q_{2}=0.419$.
Fig. 11 plots the KL divergence and energy usage for the $n$ samples for this
TD device across 10 different sets of samples. As seen Fig. 11, increasing the
number of samples lowers the KL-divergence from the desired probability
distribution. However, more samples come at the cost of increased energy
consumption. One of the resulting distributions with 2000 samples is shown in
Fig. 12.

Figure 10: Optimized weight values for each device over twenty optimization
runs using LEAP. Figure 11: KL divergence and energy usage vs. number of
samples for the given distribution with the TD device. There were 10 trials
conducted per sample. Sample sizes investigated include 10, 50, 100, 200, 500,
1000, 1500 and 2000. Figure 12: Empirical distribution using TD for 2000
samples in a single run. The optimized values generated were: $w=0.714$,
$p_{1}=0.891$, $p_{2}=0.107$, $q_{1}=0.766$, and $q_{2}=0.419$.

###  V-E Conclusions

We leveraged LEAP to generate optimal device parameters for a given
probability distribution given in equations (2). These parameters were then
plugged into device models to generate empirical distributions as shown in
Figs. 7, 9 and 12. In Fig. 10, we show the 20 different sets of weight values
that were optimized for each device type for $\omega_{1}=7500$,
$\omega_{2}=0.005$, and $\omega_{3}=0.5$. Here, we see that the weights are
customized for the device’s behavior to target the best performance in terms
of KL divergence and energy usage.

One of the challenges in optimizing for both algorithms and devices was
appropriately abstracting the device models and algorithmic constraints. The
functional models developed will also be evolved in time as new device data
and research emerges. Our framework is set up so as to accommodate such
changes.

##  VI Future Work

In this work, we leveraged LEAP for parameter optimization; however, LEAP does
not allow for network topology optimization. For future work, for both
topology and parameter optimization, we will leverage Evolutionary
Optimization for Neuromorphic Systems (EONS) [7, 8] for evolving both non-
neuromorphic and neuromorphic circuits of these devices. In this work, we
specifically explored a “Hidden Dependence” method. Future work will include
other methods to evaluate correlation and rank them based on performance
tradeoff. We will also explore the tradeoffs in circuit size, number of
devices, sampling rates, topology, and performance. We also plan to expand our
optimization technique to include system design and architecture design.

Modeling complex problems such as nuclear and high energy physics, climate
models with high precision, and novel AI techniques require simulating
probabilistic behaviors. However, doing so in conventional digital hardware
has a high energy cost. In this work, we explored drawing samples directly
from a simple, non-uniform distribution using true random number generating
(TRNG) device models. We chose a simple distribution among four outcomes that
cannot be represented by the outcome of two independent coin tosses alone.
Rather, we employed probabilistic mixing through a hidden process that can be
represented as a probabilistic circuit. Our choice of distribution was
important to demonstrate a methodology to correlate output among independent
tosses. Leveraging correlation of output is critical in direct drawing from
non-trivial distributions. We demonstrated this method using three device
models, namely MTJ-SHE, MTJ-VCMA and TD. Direct random number generation
leveraging nanoscale devices promise significant energy and latency advantages
for scientific computing applications.

Our ultimate goal is to achieve a billion random number per microsecond using
TRNG devices and computing on distributions with low latency and energy
footprint as seen in the brain [5]. The results presented in this paper are a
step in that direction. Our research intends to develop probabilistic neural
computing, which will be relevant for both AI and scientific computing
applications.

## Acknowledgment

The authors acknowledge financial support from the DOE Office of Science (ASCR
/ BES) for our Microelectronics Co-Design project COINLFIPS. Sandia National
Laboratories is a multimission laboratory managed and operated by National
Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary
of Honeywell International Inc., for the U.S. Department of Energy’s National
Nuclear Security Administration under contract DE-NA0003525.

This article has been authored by an employee of National Technology &
Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the
U.S. Department of Energy (DOE). The employee owns all right, title and
interest in and to the article and is solely responsible for its contents. The
United States Government retains and the publisher, by accepting the article
for publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this article or allow others to do so, for United States
Government purposes. The DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. SAND Number: SAND2022-16607 C.

## References

  * [1] Y. Huang and R. P. Rao, “Neurons as monte carlo samplers: Bayesian inference and learning in spiking networks,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [2] L. Buesing, J. Bill, B. Nessler, and W. Maass, “Neural dynamics as sampling: a model for stochastic computation in recurrent networks of spiking neurons,” _PLoS computational biology_ , vol. 7, no. 11, p. e1002211, 2011.
  * [3] J. D. Smith, A. J. Hill, L. E. Reeder, B. C. Franke, R. B. Lehoucq, O. Parekh, W. Severa, and J. B. Aimone, “Neuromorphic scaling advantages for energy-efficient random walk computations,” _Nature Electronics_ , vol. 5, no. 2, pp. 102–112, 2022.
  * [4] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, “Epos lhc: Test of collective hadronization with data measured at the cern large hadron collider,” _Physical Review C_ , vol. 92, no. 3, p. 034906, 2015.
  * [5] S. Misra, L. Bland, S. Cardwell, J. Incorvia, C. James, A. Kent, C. Schuman, J. Smith, and A. JB, “Probabilistic neural computing with stochastic devices,” _Advanced Materials (In Press)_ , no. 1, 2022.
  * [6] M. A. Coletti, E. O. Scott, and J. K. Bassett, “Library for evolutionary algorithms in python (leap),” in _Proceedings of the 2020 Genetic and Evolutionary Computation Conference Companion_ , 2020, pp. 1571–1579.
  * [7] C. D. Schuman, J. S. Plank, A. Disney, and J. Reynolds, “An evolutionary optimization framework for neural networks and neuromorphic architectures,” in _2016 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2016, pp. 145–154.
  * [8] C. D. Schuman, J. P. Mitchell, R. M. Patton, T. E. Potok, and J. S. Plank, “Evolutionary optimization for neuromorphic systems,” in _Proceedings of the Neuro-inspired Computational Elements Workshop_ , 2020, pp. 1–9.
  * [9] C. Mattiussi and D. Floreano, “Analog genetic encoding for the evolution of circuits and networks,” _IEEE Transactions on evolutionary computation_ , vol. 11, no. 5, pp. 596–607, 2007.
  * [10] R. S. Zebulum, M. A. Pacheco, and M. Vellasco, “Analog circuits evolution in extrinsic and intrinsic modes,” in _International Conference on Evolvable Systems_. Springer, 1998, pp. 154–165.
  * [11] J. Lehman, J. Clune, D. Misevic, C. Adami, L. Altenberg, J. Beaulieu, P. J. Bentley, S. Bernard, G. Beslon, D. M. Bryson _et al._ , “The surprising creativity of digital evolution: A collection of anecdotes from the evolutionary computation and artificial life research communities,” _Artificial life_ , vol. 26, no. 2, pp. 274–306, 2020.
  * [12] H. Schmid, C. Bessire, M. T. Bjork, A. Schenk, and H. Riel, “Silicon nanowire esaki diodes,” _Nano Letters_ , vol. 12, no. 2, pp. 699–703, 2012.
  * [13] S.-Y. Chung, R. Yu, N. Jin, S.-Y. Park, P. R. Berger, and P. E. Thompson, “Si/SiGe resonant interband tunnel diode with $f_{r0}$ 20.2 GHz and peak current density 218 $kA/cm^{2}$ for K-band mixed-signal applications,” _IEEE Electron Device Letters_ , vol. 27, no. 5, pp. 364–367, 2006.
  * [14] R. Bernardo-Gavito, I. E. Bagci, J. Roberts, J. Sexton, B. Astbury, H. Shokeir, T. McGrath, Y. J. Noori, C. S. Woodhead, M. Missous _et al._ , “Extracting random numbers from quantum tunnelling through a single diode,” _Scientific reports_ , vol. 7, no. 1, pp. 1–6, 2017.
  * [15] K. Y. Camsari, B. M. Sutton, and S. Datta, “p-bits for probabilistic spin logic,” _Applied Physics Reviews_ , vol. 6, no. 011305, pp. 1931–9401, 2019\.
  * [16] W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno, and S. Datta, “Integer factorization using stochastic magnetic tunnel junctions,” _Nature_ , vol. 573, no. 31534247, pp. 390–393, 2019.
  * [17] K. Hayakawa, S. Kanai, T. Funatsu, J. Igarashi, B. Jinnai, W. A. Borders, H. Ono, , and S. Fukami, “Nanosecond random telegraph noise in in-plane magnetic tunnel junctions,” _Physical Review Letters_ , vol. 126, no. 117202, pp. 390–393, 2021.
  * [18] C. Safranski, J. Kaiser, P. Trouilloud, P. Hashemi, G. Hu, and J. Z. Sun, “Demonstration of nanosecond operation in stochastic magnetic tunnel junctions,” _Nano letters_ , vol. 21, no. 33630604, pp. 2040–2045, 2021\.
  * [19] Z. Wang, W. Zhao, E. Deng, J.-O. Klein, and C. Chappert, “Perpendicular-anisotropy magnetic tunnel junction switched by spin-hall-assisted spin-transfer torque,” _Journal of Physics D: Applied Physics_ , vol. 48, no. 065001, pp. 065 001–1–065 001–7, 2015.
  * [20] J. Leliaert, J. Mulkers, J. D. Clercq, A. Coene, M. Dvornik, and B. V. Waeyenberge, “Adaptively time stepping the stochastic landau-lifshitz-gilbert equation at nonzero temperature: Implementation and validation in mumax3,” _AIP Advances_ , vol. 7, no. 125010, pp. 125 010–1–125 010–13, 2017.
  * [21] K. Gryszka, “From biased coin to any discrete distribution,” _Periodica Mathematica Hungarica_ , vol. 83, no. 1, pp. 71–80, 2021.