# GFN-SR: Symbolic Regression with Generative Flow Networks

Sida Li
The University of Chicago
<EMAIL_ADDRESS>
&Ioana Marinescu
Princeton University
<EMAIL_ADDRESS>
Sebastian Musslick
University of Osnabrück, Brown University
<EMAIL_ADDRESS>

###### Abstract

Symbolic regression (SR) is an area of interpretable machine learning that
aims to identify mathematical expressions, often composed of simple functions,
that best fit in a given set of covariates $X$ and response $y$. In recent
years, deep symbolic regression (DSR) has emerged as a popular method in the
field by leveraging deep reinforcement learning to solve the complicated
combinatorial search problem. In this work, we propose an alternative
framework (GFN-SR) to approach SR with deep learning. We model the
construction of an expression tree as traversing through a directed acyclic
graph (DAG) so that GFlowNet can learn a stochastic policy to generate such
trees sequentially. Enhanced with an adaptive reward baseline, our method is
capable of generating a diverse set of best-fitting expressions. Notably, we
observe that GFN-SR outperforms other SR algorithms in noisy data regimes,
owing to its ability to learn a distribution of rewards over a space of
candidate solutions.

##  1 Introduction

As a rising field that has garnered growing interest in the AI for Science
community, symbolic regression (SR) aims to solve the following problem: given
a dataset of covariates $X\in\mathbb{R}^{n\times d}$ and response
$y\in\mathbb{R}^{n}$ (where $n$ is sample size and $d$ is the dimension of
data), SR seeks to search for an interpretable expression
$f:\mathbb{R}^{d}\to\mathbb{R}$ that best explains the dataset. Specifically,
$f$ is required to be a closed-form symbolic expression that consists of
symbols from a pre-defined library $L$ containing operators like
$\exp(\cdot),\>\div$, and variables like $x_{i},1\leq i\leq d$. By learning
both the structure and parameters (i.e. constants and variables) of a
mathematical expression, SR thus explores a far richer search space and is
more flexible than traditional regression tasks (e.g. linear regression) with
fixed functional forms. On the other hand, the solution of SR is naturally
more interpretable and succinct than black-box function approximations (e.g.
non-linear neural networks), making it particularly useful in fields where
interpretability matters, such as physics or cognitive science ([15]; [12]).

Despite its appeal for automated scientific discovery, SR is widely
acknowledged to be a challenging problem ([13]; [9]). In fact, due to the
combinatorial nature of the underlying search problem, SR can be framed as an
NP-hard problem ([20]). Historically, SR has been tackled through
combinatorial optimization techniques like Genetic Programming (GP; [8];
[21]). GP maintains a population of symbolic expressions that evolves over
time to enhance the best-fit function for a given dataset. However, GP is
known to be sensitive to initialization and choice of hyper-parameters, and it
scales poorly with the data size due to its computational complexity ([4]).

In recent years, there has been a flourishing of new SR methods. Although
varying in many aspects, these methods coincidentally leverage the binary-tree
representation of symbolic expressions. Bayesian Symbolic Regression (BSR)
specifies the priors behind different operators and the likelihood of
different modifications on the expression tree ([6]). It thus approximates the
high-dimensional posterior distribution of suitable symbolic expressions for a
dataset by taking Reversible Jump MCMC tree samples. On the other hand, AI-
Feynman ([16]) relies on recursively exploiting the subtree of an expression
and solving the simplified problems.

Most notably, Deep Symbolic Regression (DSR) exploits the discrete and
sequential nature of building a valid expression tree ([13]). Under a deep
reinforcement learning framework, DSR models the construction process of an
expression tree as querying a trajectory of “tokens”
$\tau=[\tau_{1},...,\tau_{T}]$ (elements from the library $L$) one at a time
from the categorical distribution
$p(\tau_{i}|\>\theta,\tau_{1},...,\tau_{i-1})$, which derives from a neural
network policy parametrized by $\theta$ (such as an RNN). The completed
expression tree will then be evaluated based on its fit to the dataset (e.g.
normalized RMSE), and the evaluation will be transformed into a reward
function $R(\tau)$ (e.g. $1/(1+\textit{NRMSE}(\tau))$) to “reward” the
trajectory. Ultimately, the objective is to optimize $\theta$ to maximize the
objective: 111In fact, the DSR paper adopted a risk-seeking version of this
objective $J_{risk}(\theta)=\mathbb{E}_{\tau\sim
p(\tau|\theta)}[R(\tau)|R(\tau)>R_{\epsilon}]$ ($\epsilon$ is a risk
threshold) and demonstrated its superiority over vanilla policy gradient (VPG)
that uses the objective above. Nevertheless, both objectives seek to maximize
an expectation and suit the needs of this paper.

$J(\theta)=\mathbb{E}_{\tau\sim p(\tau|\theta)}[R(\tau)]$ (1)

i.e. the expected reward along the trajectory $\tau$, which translates to the
expected fitness to the dataset for an expression generated by DSR.

Despite proving their effectiveness on SR benchmarks, DSR and its variants
suffer from an inherent limitation within the RL setup: maximization of
expected reward is often achieved through exploiting a single high-reward
sequence of actions ([14]). The issue becomes particularly critical under
noisy settings, where equations with symbolic differences might produce
datasets that appear deceptively similar. Therefore, the ability to generate
diverse high-reward candidate equations—instead of clinging onto a single
highest-reward solution—is a desirable property in such scenarios ([1]).

In this paper, we present GFN-SR, a deep-learning SR method based on
Generative Flow Networks (GFlowNets; GFN; [2]). GFN-SR sequentially samples
tokens from a stochastic policy (represented by a deep neural net) to
construct expression trees, similar to DSR. However, by casting SR into a
GFlowNets problem instead of a return-maximization RL one, GFN-SR aims to
align the distribution over the generated expressions proportional to their
fit to the dataset, measured by a reward function $R(\cdot)$. Therefore, GFN-
SR is able to not only identify a high-reward expression when one exists but
also generate a set of diverse candidates when the reward distribution over
the space of expressions is multi-modal.

To summarize our contributions: (1) we propose a novel deep-learning SR method
by leveraging GFlowNets as sequential samplers, (2) we come up with an
adaptive reward function that balances exploration with exploitation, (3) we
report experiments evaluating our method’s performance in settings with and
without noise, showcasing its superior performance over other methods when the
data are noisy.

##  2 Methods

###  2.1 SR under GFlowNets framework

Figure 1: A part of the state space $\mathcal{S}$ and transitions between
expression trees when we limit the maximum depth to two (see discussion of
depth constraint below). The source state $s_{0}$ is an empty tree on the
left, and terminal states are complete expression trees within orange plates.
Any node in grey indicates that it has not been constructed yet. The token
library $L$ includes, but is not limited to, $sin$, $\times$, and variables
$x,y$.

In this section, we formalize the Symbolic Regression problem within the
GFlowNets framework ([2]). We denote a complete expression tree as a binary
tree representation of mathematical expressions with non-leaf nodes being
operators or functions (in this paper, we use “operators” specifically for
binary operators and “functions” for unary ones), and leaf nodes being
variables and constants. Each token in the expression tree is part of a pre-
defined library $L$. A complete expression tree $s$ can be reproduced from
scratch: starting with an empty tree $s_{0}$, one can take the pre-order
traversal of $s$ and sequentially add nodes to the tree. We can define
$\mathcal{S}$ to be the overall set of complete expression trees and all
intermediate trees in the above sequential construction process. Similarly, we
denote $\mathcal{A}$ as the set of directed edges that represent all possible
transitions between elements in $\mathcal{S}$ towards constructing some
complete expression tree $s\in\mathcal{S}$. Trivially, since a transition can
only happen when we add a new node to the previous intermediate tree, the
combination $\\{\mathcal{S},\mathcal{A}\\}$ creates a directed acyclic graph
(DAG) for GFlowNets to be trained on (in fact, $\\{\mathcal{S},\mathcal{A}\\}$
is a tree graph due to the pre-order traversal specification). In this DAG,
$s_{0}$ is the source and any complete expression tree is a sink (terminal
state). Figure 1 illustrates a small part of the state space $\mathcal{S}$ of
expression trees and the possible transitions from an empty tree to a valid
expression.

For the given DAG $\\{\mathcal{S},\mathcal{A}\\}$ as well as a reward function
$R(\cdot)$ acting on the terminal states in $\mathcal{S}$, GFN is a generative
model with the objective to sequentially generate a complete expression tree
$s\in\mathcal{S}$ with probability ([1])

$\pi(s)\approx\frac{R(s)}{Z}=\frac{R(s)}{\sum_{s^{\prime}\in\mathcal{S}^{c}}R(s^{\prime})}$
(2)

where $Z$ is a normalizing constant summing up all terminal rewards, and
$\mathcal{S}^{c}\subset\mathcal{S}$ is the set of all complete expression
trees. Due to the size of $\mathcal{S}^{c}$ being often intractable, GFN
treats $Z$ as a learnable parameter. To achieve the objective in Equation 2,
GFN also needs to learn a stochastic forward-sampling policy
$P_{F}(\cdot|\cdot)$, which is represented by a deep neural network with
parameter $\theta$. For a given state $s\in\mathcal{S}$, $P_{F}(\cdot|s)$
describes a categorical distribution with support
$\\{\tilde{s}\in\mathcal{S}:s\to\tilde{s}\in\mathcal{A}\\}$. In our context,
$P_{F}(\cdot|s)$ boils down to a categorical distribution over all possible
tokens from $L$ and governs the next construction step from $s$.

To optimize parameters $\theta$ for the stochastic policy in GFN-SR, we adopt
the trajectory balance (TB) loss proposed by [11]. Let $\tau=(s_{0}\to
s_{1}\to\cdots\to s_{n}=s)$ be the trajectory of generating a complete
expression tree $s$ from scratch, we have:

$\mathcal{L}_{TB}(\tau)=\Bigg{(}\log\frac{Z_{\theta}\prod_{i=1}^{n}P_{F}(s_{t}|s_{t-1};\theta)}{R(s)\prod_{i=1}^{n}P_{B}(s_{t-1}|s_{t})}\Bigg{)}^{2}$
(3)

where we condition both $Z$ and $P_{F}$ on $\theta$ to highlight their
dependencies. The $P_{B}(s_{t-1}|s_{t})$ term appearing in the denominator
refers to a backward-sampling policy, which estimates the distribution of
predecessors for a given state in the construction and is usually jointly
learnt with $P_{F}$. Fortunately, the tree structure of
$\\{\mathcal{S},\mathcal{A}\\}$ in our problem setup gives each non-empty
expression tree in $\mathcal{S}$ a unique predecessor (we can simply identify
the most recently inserted node by reverting the pre-order and remove this
node to obtain the predecessor). The TB loss for GFN-SR thus simplifies to:

$\mathcal{L}_{TB}^{*}(\tau)=\Bigg{(}\log\frac{Z_{\theta}\prod_{i=1}^{n}P_{F}(s_{t}|s_{t-1};\theta)}{R(s)}\Bigg{)}^{2}$
(4)

[11] has shown that minimizing the TB loss is consistent with achieving the
desirable objective in Equation 2. Finally, the stochastic policy is optimized
via mini-batch gradient descent on $\theta$.

Figure 2: Illustration of the sequential expression tree construction. Through
stochastic policy emits a raw categorical distribution, which is adjusted and
re-normalized by the in-situ constraints (probabilities being zeroed out are
labelled in red color). Finally, tokens for the next tree nodes are sampled
and added to the growing expression tree following pre-order. This figure also
corresponds to Figure 1, which views the same process but from a state space
perspective.

###  2.2 LSTM as the policy network

A key design choice in GFlowNets is the deep neural network behind the
stochastic (forward) policy. As our goal is to sequentially predict the next
token (symbol in our library $L$), we take inspiration from recurrent neural
network (RNN) models that have been widely adopted by deep SR methods for
similar sequence prediction tasks ([13], [10]). In particular, GFN-SR
leverages Long Short-Term Memory (LSTM) networks ([5]). Through the inclusion
of gating mechanisms (input, forget, and output gates) for the memory cell
state, LSTMs can learn longer-term dependencies in the input sequences ([18]).
This facet of LSTMs aligns with our need to construct subsequent segments of
expression trees based on prior, and potentially distant, tokens. In our
implementation, we feed the one-hot representation of parent and sibling
tokens of the next token to be generated into our LSTM policy network 222due
to the pre-order, the parent token exists for any node but the root, and
sibling exist for all nodes as the right child. We use a placeholder empty
token if either does not exist.. The LSTM then emits a probability vector of
length equal to the size of library $L$ that represents a categorical
distribution over the next token. Details on our LSTM architecture and hyper-
parameters are listed in the Appendix.

###  2.3  Applying constraints to reduce computational complexity

Before a token is sampled out based on the probabilities emitted by the policy
network, we apply in-situ constraints on the space of possible tokens by
zeroing out certain entries in the probability vector and reweighing it. These
constraints effectively reduce the computational complexity and accelerate the
search for high-quality candidates. They are rules based on realistic
assumptions that fall into the following categories:

Composition constraints We disallow certain compositions of operators and
functions since they rarely exist in real equations ([13]). For instance, a
trigonometric function can not be nested in another one so expressions like
$\sin(\cos(x)+\sin(x))$ is ruled out. We also prohibit functions that are
inverse with each other to be composed directly (e.g. $(\sqrt{x})^{2}$).

Depth constraints We constrain how long an expression can be by imposing
constraints on the depth (or height) of its corresponding expression tree.
When the next node of a tree under construction is at maximum depth, we zero
out all probabilities except for tokens representing variables and the
constant symbol. This constraint is especially useful since we often process a
batch of expressions together and the longest expression dictates the
computation speed through the policy network.

Constant constraints Any function (i.e. unary operator) cannot take constant
as its argument since that would simply represent another constant; any binary
operator can only take constant as its right child, which (1) reduces
expression redundancies due to the symmetry of tree structure and (2) prevents
both arguments of such operator to be constants. The constant tokens are
specially handled in GFN-SR as they are not assigned specific values during
construction, and we instead optimize them in a separate process after full
expressions are built (details in Appendix C). These constraints thus save a
considerable amount of time fitting constants.

Figure 2 uses a concrete example to illustrate the two procedures (2.2, 2.3)
above in sequentially generating an expression tree.

###  2.4 Reward baseline for GFlowNets

In this section, we discuss the role of reward function $R$ and a novel design
of reward baseline for GFlowNets (the terminology “baseline” is inspired by
[22]).

Based on Equation 2, a properly trained GFlowNet learns to generate a terminal
state $s\in\mathcal{S}$ with probability $\pi(s)\propto R(s)$, thus the reward
function should reflect our satisfaction on choosing $s$ as a solution to our
problem. For the SR task, a reward function assigns a non-negative scalar to a
terminal state (expression tree) $s\in\mathcal{S}$ to measure how well $s$
fits in the given dataset $(X,y)$. Usually, the reward function is some
transformation on a loss function (e.g. RMSE, MAE) that quantifies such
fitness. We assume $R(s)=1/(1+RMSE(s))$ as the “vanilla” reward function in
the following discussion.

Generating a diverse set of high-quality candidate expressions requires a good
balance between exploration and exploitation. GFlowNets are good at
exploration by design since their training objective encourages generating
expressions with probabilities proportional to the rewards. However, if the
landscape of $R(s)$ is multi-modal yet the modes are not peaked enough,
GFlowNets would find it difficult to concentrate on regions with high rewards.
To address this issue, we can adjust $R$ such that low-reward expressions
receive even lower rewards, while high-performing solutions are awarded
greater rewards to emphasize their superiority. As a result, the distribution
of $R$ would concentrate around the modes. Denote $B>0$ as a baseline that
does not depend on $s$, and consider:

$R_{B}(s)=R(s)\bigg{(}1+\gamma\Big{(}\frac{R(s)}{B}-1\Big{)}\bigg{)}$ (5)

where $\gamma\in[0,1]$ is a scaling hyper-parameter that controls the degree
of adjustment. We can easily verify that $R_{B}(s)$ is non-negative and
satisfies the aforementioned properties. In the actual implementation, $B$ is
initialized to be the average vanilla reward of expressions in the first
batch. As training proceeds, we gradually increase $B$ toward the average
vanilla reward for the top-performing expressions seen so far. The pseudo-
codes for calculating the baseline and adjusted reward are included in
Appendix A. Source code for GFN-SR will be available at
https://github.com/listar2000/gfn-sr.

##  3 Experiments

###  3.1 Noiseless Nguyen benchmark

We first experiment GFN-SR on the Nguyen benchmark problems ([17]) to evaluate
its general capacity in fitting the given dataset (in terms of RMSE) when
there is no noise. The Nguyen benchmark includes 12 distinct expressions of
various complexity (see below) and has been widely used to evaluate SR
methods.

$\displaystyle h_{1}$ $\displaystyle=x^{3}+x^{2}+x$ $\displaystyle h_{2}$
$\displaystyle=x^{4}+x^{3}+x^{2}+x$ $\displaystyle h_{3}$
$\displaystyle=x^{5}+x^{4}+x^{3}+x^{2}+x$ $\displaystyle h_{4}$
$\displaystyle=x^{6}+x^{5}+x^{4}+x^{3}+x^{2}+x$ $\displaystyle h_{5}$
$\displaystyle=\sin(x^{2})\cos(x)-1$ $\displaystyle h_{6}$
$\displaystyle=\sin(x)+\sin(x+x^{2})$ $\displaystyle h_{7}$
$\displaystyle=\log(x+1)+\log(x^{2}+1)$ $\displaystyle h_{8}$
$\displaystyle=\sqrt{x}$ $\displaystyle h_{9}$
$\displaystyle=\sin(x)+\sin(y^{2})$ $\displaystyle h_{10}$
$\displaystyle=2\sin(x)\cos(y)$ $\displaystyle h_{11}$ $\displaystyle=x^{y}$
$\displaystyle h_{12}$ $\displaystyle=x^{4}-x^{3}+\frac{y^{2}}{2}-y$

For comparison, we also run the same benchmarks on three other SR methods: GP,
DSR, and BSR. Due Each of these methods have their unique hyper-parameters and
architectures, and we tune them deliberately so that the total number of
expression evaluations in one run (for one benchmark function) is around 25
million. We perform 20 repeating trials (but with different seeds) of the
Nguyen benchmarks for each method and calculate the mean & standard error of
the RMSEs (for detailed experimental setups for each method, see Appendix C).

Table 1: RMSE comparison of SR methods on Nguyen benchmarks
(standard error in parenthesis)) Benchmark | Dataset | GFN-SR | DSR | GP | BSR
---|---|---|---|---|---
Nguyen-1 | U(20, -1, 1) | $0(0)$ | $0(0)$ | $0(0)$ | $0(0)$
Nguyen-2 | U(20, -1, 1) | $0(0)$ | $0(0)$ | $0(0)$ | .002(.001)
Nguyen-3 | U(20, -1, 1) | $0(0)$ | $0(0)$ | $0(0)$ | .01(.004)
Nguyen-4 | U(20, -1, 1) | $.02(.002)$ | 0(0) | .05(.005) | .01(.003)
Nguyen-5 | U(20, -5, 5) | .15(.021) | .35(.028) | .51(.023) | .19(.066)
Nguyen-6 | U(20, -5, 5) | .25(.082) | .012(.005) | .90(.435) | .47(.158)
Nguyen-7 | U(20, 0, 2) | .128(.035) | .014(.0002) | .044(.002) | 6.2e-4 (2.9e-4)
Nguyen-8 | U(20, 0, 4) | 0(0) | 0(0) | .091(0) | 9.6e-4 (4.6e-4)
Nguyen-9 | U(20, $[0,1]^{2}$) | .023(.015) | 0(0) | .048(.011) | .001(.001)
Nguyen-10 | U(20, $[0,1]^{2}$) | .053(.014) | .013(.019) | .044(.006) | .025(.014)
Nguyen-11 | U(20, $[0,1]^{2}$) | .025(.012) | .017(.021) | .069(.012) | .054(.005)
Nguyen-12 | U(20, $[0,1]^{2}$) | .037(.014) | .040(.004) | .080(.021) | .082(.003)

In general, results in 1 demonstrate the superiority of deep SR methods in
terms of producing best-fitting solutions to a given dataset without noise.
While DSR takes the lead for many benchmarks, we are thrilled to observe that
GFN-SR, which is not trained under a return-maximization objective, still
exhibits competitive performance and even outperforms other methods in
Nguyen-5 and Nguyen-12. These results demonstrate the potential and efficacy
of our new framework leveraging GFlowNets in symbolic regression tasks.

###  3.2 Noisy synthetic dataset

We set up a synthetic benchmark that is specially designed to evaluate an SR
method’s ability to recover exact symbolic expressions under noise. The
benchmark consists of six equations that are grouped into three pairs (we
denote equations in the $i$th pair as $f_{i}$ and $g_{i}$):

Pair 1: $\displaystyle\quad f_{1}=x^{3}+x^{2}+x\quad
g_{1}=\sin(x)(\sqrt{x}+\exp{x})\quad x\in[0,1]$ Pair 2: $\displaystyle\quad
f_{2}=\log(x+1)\quad\>\>g_{2}=0.5\cdot x\quad x\in[-0.1,0.1]$ Pair 3:
$\displaystyle\quad f_{3}=x^{2}-4x+3\quad g_{3}=4\cdot(1-\sqrt{x})\quad
x\in[0,1]$

These equations are seemingly simple and solvable in noiseless settings (in
fact $f_{1}$ is the first equation in the Nguyen benchmark). However, within
the domains we specify, equations in each group are symbolically different yet
very close to each other (see Appendix B for visualization). Therefore, if we
further add noises to the dataset, it becomes challenging for an SR method to
recover the right equation, especially if the method lacks the ability to
explore a wide region of solutions and easily sticks around a single high-
reward mode.

For each equation, the dataset consists of 40 points uniformly sampled from
its specified domain (20 for training and 20 for testing) and responses $y$
under 10% noise level (i.e. noises $\epsilon\sim 0.1\cdot N(0,Var(y))$. For
each method (GFN-SR, DSR, and BSR), we repeat the experiment for 20 trials
with different seeds and make sure that the number of equation evaluations in
each trial is around 25 million for all methods. The recovery rate is defined
as the percentage of trails (among all 20 runs) in which exact symbolic
equivalence is achieved by the top predicted equations of a method (for DSR
and BSR, this is achieved if any equation on the Pareto frontier matches the
ground-truth; for GFN-SR, we compare the ground-truth with the most frequently
occurring equations in the samples post-training, which is a stricter
condition).

Table 2: Recovery rate of the noisy benchmark under 10% noise | Group 1 | Group 2 | Group 3
---|---|---|---
Equation | $f_{1}$ | $g_{1}$ | $f_{2}$ | $g_{2}$ | $f_{3}$ | $g_{3}$
Dataset | $U(40,0,1)$ | $U(40,-0.1,0.1)$ | $U(40,0,1)$
BSR | 30% | 5% | 0% | 90% | 5% | 75%
DSR | 75% | 5% | 5% | 40% | 5% | 90%
GFN-SR | 100% | 15% | 15% | 95% | 95% | 85%

Results in Table 2 show how SR methods are sensitive to the presence of noise,
as the recovery rates are consistently low even for basic expressions like
$f_{2}$. Bayesian Symbolic Regression (BSR) represents equations as a linear
mixture of individual trees ([6]), thus it is not good at recovering equations
exactly. DSR’s return maximization objective, as mentioned earlier, sometimes
misleads it to focus on the wrong but high-reward mode (i.e. the other
expression in the group) while overlooking the correct one. Overall, GFN-SR is
capable of generating a diverse set of candidate expressions and has
outperformed the two other methods in recovering every equation except
$g_{3}$.

##  4 Conclusion and future work

In this paper, we present GFN-SR—a novel symbolic regression method based on
GFlowNets. The GFlowNets objective enables our method to explore and generate
diverse candidate expressions, while the introduction of a reward baseline
enables us to exploit high-reward regions only in the reward landscape. These
exciting features are demonstrated through numerical experiments, where GFN-SR
shows competitive performance in noiseless benchmarks and a superior ability
to recover expressions under noise.

However, the proven success of GFN-SR is limited to experiments on the Nguyen
benchmarks and our synthetic dataset. To further verify the effectiveness and
understand the numerical properties of our approach, an immediate next step is
to perform experiments on larger datasets, such as SRBench ([9]) and Feynman
SR database ([16]). We would also want to provide noisy data from real-world
problems to confirm GFN-SR’s exciting potential in recovering expressions from
noise.

Another limitation of our current method is the enforcement of pre-order
traversal in the tree construction process. Despite some aforementioned
benefits (e.g. simplification of backward policy $P_{B}$), this restriction
severely limits the number of possible transitions in the state space. If we
allow growing any part of an incomplete tree, $\\{\mathcal{S},\mathcal{A}\\}$
will still be a valid DAG that we can train GFlowNets on. This more flexible
approach could potentially enhance the exploration of the expression space,
though it may require additional computational time.

Finally, our presented method is versatile, and many of its components admit
new ideas. For instance, Transformers have recently been the prevailing deep
learning model for sequence generation tasks in NLP ([19]; [3]). They have the
potential to replace RNNs as policy networks and provide a better
understanding of global structures within symbolic expressions. Also, [10]
proposed a unified framework to combine deep-learning-based SR methods with
other techniques, such as GP and polynomial regression. It is worth
investigating to see how we can further improve GFN-SR under this framework.

## Acknowledgement

We thank Younes Strittmatter, Chad Williams, Ben Andrew, and other members of
the Autonomous Empirical Research Group for their feedback and insights in
this work. We appreciate the computation resources provided by the Center for
Computation & Visualization (CCV) at Brown University. Sebastian Musslick was
supported by Schmidt Science Fellows, in partnership with the Rhodes Trust,
and supported by the Carney BRAINSTORM program at Brown University.

## References

  * [1] Emmanuel Bengio et al.  “Flow network based generative models for non-iterative diverse candidate generation”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 27381–27394
  * [2] Yoshua Bengio et al.  “Gflownet foundations”  In _arXiv preprint arXiv:2111.09266_ , 2021
  * [3] Anthony Gillioz, Jacky Casas, Elena Mugellini and Omar Abou Khaled  “Overview of the Transformer-based Models for NLP Tasks”  In _2020 15th Conference on Computer Science and Information Systems (FedCSIS)_ , 2020, pp. 179–183  IEEE
  * [4] Ahmad B Hassanat et al.  “An improved genetic algorithm with a new initialization mechanism based on regression techniques”  In _Information_ 9.7 MDPI, 2018, pp. 167
  * [5] Sepp Hochreiter and Jürgen Schmidhuber  “Long short-term memory”  In _Neural computation_ 9.8 MIT press, 1997, pp. 1735–1780
  * [6] Ying Jin et al.  “Bayesian Symbolic Regression”, 2020  arXiv:1910.08892 [stat.ME]
  * [7] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [8] John R Koza  “Hierarchical genetic algorithms operating on populations of computer programs.”  In _IJCAI_ 89, 1989, pp. 768–774
  * [9] William La Cava et al.  “Contemporary symbolic regression methods and their relative performance”  In _arXiv preprint arXiv:2107.14351_ , 2021
  * [10] Mikel Landajuela et al.  “A unified framework for deep symbolic regression”  In _Advances in Neural Information Processing Systems_ 35, 2022, pp. 33985–33998
  * [11] Nikolay Malkin et al.  “Trajectory balance: Improved credit assignment in gflownets”  In _Advances in Neural Information Processing Systems_ 35, 2022, pp. 5955–5967
  * [12] Sebastian Musslick  “Recovering Quantitative Models of Human Information Processing with Differentiable Architecture Search.”  In _Proceedings of the Annual Meeting of the Cognitive Science Society_ 43 eScholarship, 2021  URL: https://escholarship.org/uc/item/9wd571ts
  * [13] Brenden K Petersen et al.  “Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients”  In _arXiv preprint arXiv:1912.04871_ , 2019
  * [14] Richard S Sutton, David McAllester, Satinder Singh and Yishay Mansour  “Policy gradient methods for reinforcement learning with function approximation”  In _Advances in neural information processing systems_ 12, 1999
  * [15] Wassim Tenachi, Rodrigo Ibata and Foivos I Diakogiannis  “Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws”  In _arXiv preprint arXiv:2303.03192_ , 2023
  * [16] Silviu-Marian Udrescu et al.  “AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 4860–4871
  * [17] Nguyen Quang Uy et al.  “Semantically-based crossover in genetic programming: application to real-valued symbolic regression”  In _Genetic Programming and Evolvable Machines_ 12 Springer, 2011, pp. 91–119
  * [18] Greg Van Houdt, Carlos Mosquera and Gonzalo Nápoles  “A review on the long short-term memory model”  In _Artificial Intelligence Review_ 53 Springer, 2020, pp. 5929–5955
  * [19] Ashish Vaswani et al.  “Attention is all you need”  In _Advances in neural information processing systems_ 30, 2017
  * [20] Marco Virgolin and Solon P Pissis  “Symbolic regression is np-hard”  In _arXiv preprint arXiv:2207.01018_ , 2022
  * [21] Ekaterina Yurievna Vladislavleva  “Model-based problem solving through symbolic regression via pareto genetic programming”  CentER, Tilburg University Tilburg, 2008
  * [22] Lex Weaver and Nigel Tao  “The optimal reward baseline for gradient-based reinforcement learning”  In _arXiv preprint arXiv:1301.2315_ , 2013

## Appendix A: Pseudocodes

Algorithm 1 Calculate reward and update baseline

1:a batch of expressions $\\{s_{1},\dots,s_{N}\\}$

2:vanilla reward function $R$

3:current baseline $B$

4:a priority queue of best vanilla rewards $Q$

5:

6:Hyper-parameters: scale parameter $\gamma$, moving average weight $\alpha$

7:$\triangleright$ calculate the vanilla and baseline-aware reward

8:for $i$ in $1$ to $N$ do

9: $R_{i}\leftarrow R(s_{i})$

10: $R^{*}_{i}\leftarrow R(s_{i})\times(1+\gamma\times(\frac{R(s_{i})}{B}-1))$

11:end for

12:$\triangleright$ update the priority queue with rewards sampled in this
batch

13:$Q\leftarrow\text{updatePQ}(Q,\\{R_{1},\dots,R_{N}\\})$

14:$\triangleright$ take the average of rewards in $Q$

15:$B^{\prime}\leftarrow\text{calcAvgReward}(Q)$

16:$B^{*}\leftarrow\alpha\times B+(1-\alpha)\times B^{\prime}$

17:

18:rewards $\\{R^{*}_{1},\dots,R^{*}_{N}\\}$, new baseline $B^{*}$

## Appendix B: Dataset

Below is the plot of functions $f_{1}(x)=x^{3}+x^{2}+x$ and
$g_{1}(x)=\sin(x)(\sqrt{x}+\exp(x))$ mentioned in Section 3.2, as well as the
noisy observations for 20 sample points at noise level 10% (so they are
similar to the training/testing dataset).

## Appendix C: Experiment Details

#### Hyper-parameter tuning

For GFN-SR, we consider three major hyper-parameters: (1) batch size
$b\in[250,500,1000]$, learning rate $\alpha\in$ [1e-4, 2e-4, 5e-4, 1e-3], and
the scaling parameter for reward baseline $\gamma\in[0.1,0.2,\cdots,0.9]$. We
tuned them on a simplified Nguyen benchmark and the final hyper-parameters
used are $b=250,\alpha=\text{5e-4},\gamma=0.9$.

#### LSTM architecture

The LSTM policy network we use throughout the paper has 2 recurrent layers and
250 hidden states for each layer. We train the LSTM using Adam optimizer with
the learning rate tuned above ([7]). Before we feed the parent and sibling
tokens into the LSTM, they would go through a one-hot encoding.

#### Token library

Throughout all experiments for all methods, the library $L$ includes functions
$\sqrt{\cdot},\log,\exp,\sin,\cos,(\cdot)^{2}$, operators $\times,\div,+,-$,
variables $x_{1},\cdots,x_{d}$ for input dimension $d$, and the constant
token. Certain SR methods require specifying the range of constants, which we
specify as the maximum and minimum constant values in our benchmarks.

#### Constant optimization

In GFN-SR, a complete expression with at least one constant token cannot be
evaluated by reward function $R$ directly. This is because these constants are
not assigned any specific value during construction; instead, a constant
optimization mechanism is applied. We have a dictionary and a counter. If an
expression with constants has not been discovered before, or if its counter is
at zero, we use BFGS (a numerical optimization algorithm) to find values of
the constants that maximize $R$. We then cache its reward into the library and
set its counter to some number that determines how frequently should we redo
the optimization step. Otherwise, we simply read the cached reward for the
expression from the dictionary and reduce the counter by one.

#### DSR

In our experiments, we use the deep symbolic optimization (DSO) Python package
to run DSR method. Specifically, DSO improves the vanilla DSR method in [13]
by providing certain post-processing. We use the same hyper-parameter
configuration that [13] used in their paper.

#### BSR

[6] have provided an official implementation of BSR; however, it is not
continually maintained and certain results from the paper cannot be reproduced
properly through this repository. We, therefore, use a third-party’s
implementation for the experiments. For BSR, we use a uniform prior for tokens
in the library and a tree num $K=3$.

#### GP

We use GPLearn, an open-source Python genetic programming toolkit, for bench-
marking SR problems with GP. The specific parameters that I used are
SymbolicRegressor(population_size=5000, const_range=(-1, 2), generations=40,
stopping_criteria=0.001, p_crossover=0.7, p_subtree_mutation=0.1,
p_hoist_mutation=0.05, p_point_mutation=0.1, max_samples=1, verbose=1,
init_method="full", parsimony_coefficient=0.01, random_state=0) and the
documentations is available here.

#### Computation resource

We run the above experiments using NVIDIA’s Quadro RTX 6000 GPU and Intel’s
Xeon Gold 6242 CPU.


# Orthogonal Spin Current Injected Magnetic Tunnel Junction for Convolutional
Neural Networks

Venkatesh Vadde1, Bhaskaran Muralidharan1,1, Abhishek Sharma2,2
<EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Electrical
Engineering, Indian Institute of Technology Bombay, Powai, Mumbai-400076,
India  2Department of Electrical Engineering, Indian Institute of Technology
Ropar, Rupnagar, Punjab-140001, India

###### Abstract

We propose that a spin Hall effect driven magnetic tunnel junction device can
be engineered to provide a continuous change in the resistance across it when
injected with orthogonal spin currents. Using this concept, we develop a
hybrid device-circuit simulation platform to design a network that realizes
multiple functionalities of a convolutional neural network. At the atomistic
level, we use the Keldysh non-equilibrium Green’s function technique that is
coupled self-consistently with the stochastic Landau-Lifshitz-Gilbert-
Slonczewski equations, which in turn is coupled with the HSPICE circuit
simulator. We demonstrate the simultaneous functionality of the proposed
network to evaluate the rectified linear unit and max-pooling functionalities.
We present a detailed power and error analysis of the designed network against
the thermal stability factor of the free ferromagnets. Our results show that
there exists a non-trivial power-error trade-off in the proposed network,
which enables an energy-efficient network design based on unstable free
ferromagnets with reliable outputs. The static power for the proposed ReLU
circuit is $0.56\mu W$ and whereas the energy cost of a nine-input rectified
linear unit-max-pooling network with an unstable free ferromagnet($\Delta=15$)
is $3.4pJ$ in the worst-case scenario. We also rationalize the magnetization
stability of the proposed device by analyzing the vanishing torque gradient
points.

##  I Introduction

Neuromorphic computing takes inspiration from biological brains to perform
highly complex problems while consuming remarkably low energy [1, 2]. The
brain performs in-memory computation and uses many low-precision calculations
in parallel to perform a task. Meanwhile, modern computers are primarily based
on the von Neumann architecture, which separates the computation and memory
units and uses high-precision calculations [3]. Convolutional neural networks
(CNN) are a class of artificial neural networks (ANNs) [4] which produces
excellent performance in machine learning problems dealing with image data
[5], computer vision [6, 7], and natural language processing [8].
A CNN has two stages: the first one is the feature extraction stage, and the
second is the fully connected stage. Feature extraction is achieved by
employing various layers, such as the convolutional layer, the activation
function layer, and the pooling layer. These layers enable the CNN to exploit
the spatial features of an image, introduce non-linearity [9], and suppress
noise[10].
The pooling layer is responsible for reducing the feature’s size, which helps
reduce the parameters and computation in the network. There are two types of
pooling: max pooling and average pooling. Max pooling calculates the maximum
value from the portion of the image feature, while the average pooling finds
the average value of the part of the image [11, 4]. The max-pooling also
performs noise suppression by discarding noisy activations and is more often
used in the network compared to the average pooling[10].
The activation function introduces non-linearity to the network. Non-linearity
enables the network to learn complex structures in the data and differentiates
between outputs [9]. Traditionally, sigmoid and tanh activation functions have
been widely utilized. But, sigmoid and tanh functions saturate when the input
is very high or low and are only sensitive to changes around their mid-points.
After saturation, the network won’t be able to learn well [10, 12]. The
sigmoid and tanh functions also face the vanishing gradient problem, where the
gradient information used to learn networks vanishes for deep networks.
Without any helpful gradient information, deep networks won’t be able to learn
effectively [10]. The vanishing gradient and saturation problems faced by tanh
and sigmoid functions can be overcome by the rectified linear unit (ReLU)
activation function[10].

Figure 1: Design schematics. (a) the MTJ device on top of a SHE layer. Charge
currents are injected along $\hat{x}$ and $\hat{y}$ directions, which produce
a spin current along the $\hat{z}$-direction with spin-polarizations along
$\hat{y}$ and -$\hat{x}$ directions respectively. (b) Circuit diagram for ReLU
output. The MTJ and the resistor $R_{1}$ form a voltage divider whose output
is fed to the CMOS inverter to produce the ReLU output. (c) Circuit diagram
for local max-pooling function comprising two ReLU circuits. The inverted
outputs of these ReLU circuits are fed as inputs to all other circuits, thus
enabling the competition required for the max-pooling functionality.

The hardware implementation of various layers of the feature extraction stage
of a CNN has been explored by a few studies [13, 14, 15], with the ultimate
aim of developing neuromorphic computing [2, 16]. These works are based on
CMOS circuit realizations of only the activation function [13, 14, 15] and
lack the concatenation ability to perform simultaneous max-pooling.
Furthermore, the in-memory computation using the existing CMOS technology is
severely limited by area and energy requirements[3, 17]. Spintronics, on the
other hand, provides a wide range of devices and physical effects for in-
memory computing that suits the hardware realization of neuromorphic computing
and enables the paradigm of “let the physics do the computing” [18]. Although
the spintronic implementation of activation functions is studied [19, 20, 21,
22] but the implementation of max-pooling is still elusive, and moreover, some
of these works require an external magnetic field to obtain the activation
function. We propose a technologically relevant spintronic implementation of
ReLU and max-pooling network using orthogonal spin current injection in MTJ
device and show that unstable ferromagnets can be used for energy-efficient
design employing an atom-to-circuit approach that weaves quantum transport
modeling with CMOS circuit design.

Most of the works centered around hardware realization of neuromorphic
computing are implementing a part of the network on hardware [23, 24, 25, 26,
27, 28, 29] (like synapses based on domain wall, MTJ, nano-oscillators) and
employ software-implementation (for activation function and pooling layers) to
complete the network. The lack of spintronic hardware implementation of the
activation function and pooling layers is a ramification of the absence of
continuous change in the resistance of the MTJ device using spin current. Our
work fills this critical gap by providing hardware for activation and max-
pooling functions, thus enabling a fully-hardware implementation of
neuromorphic networks.

The paper is organized as follows. In Sec. II, we describe the basics of the
ReLU and the max-pooling functionalities. In Sec.III, we discuss in detail our
developed hybrid simulation platform based on Keldysh non-equilibrium Green’s
function (NEGF) formalism coupled with the Landau-Lifshitz-Gilbert-Slonczewski
(LLGS) equations for quantum transport description at the device level, which
is then interfaced with the HSPICE circuit simulator to capture the spin and
charge current interplay along with CMOS devices. In Sec.IV we present the
performance of the ReLU and max-pooling circuit in terms of power consumption
and noise analysis. We show that there exists a non-trivial power-error trade-
off in the proposed ReLU-max-pooling network, which enables an energy-
efficient network design based on unstable free ferromagnets with reliable
outputs. We conclude in Sec. V.

##  II Design

###  II-A ReLU circuit

The ReLU function emulation requires a device with linear characteristics. We
propose a spin Hall effect (SHE)-driven MTJ device with orthogonal currents,
as shown in Fig. 1(a), to generate a linear and continuous rotation in
magnetization of the free FM layer. In this work, we have used typical CoFeB-
based FMs for the fixed and free layers with the equilibrium magnetizations
along the $\hat{y}$-direction and the $\hat{z}$-directions, respectively. The
fixed FM has in-plane magnetic anisotropy, and the free FM has a perpendicular
magnetic anisotropy (PMA) [30] in order to achieve the respective equilibrium
magnetization directions. The SHE layer generates two spin currents with
orthogonal polarizations corresponding to the direction of input charge
currents. These spin currents interact with the free FM layer to rotate the
magnetization to the in-plane direction. The continuous rotation in
magnetization can be translated to the electrical signal using the TMR of the
MTJ device. In order to achieve the ReLU functionality, the SHE-driven MTJ is
connected to a resistor to form a voltage divider, as shown in Fig. 1(b). The
voltage divider drives a CMOS inverter pair that operates in the linear region
to produce the ReLU output.
Continuous-linear rotation in the magnetization of the free FM layer can also
be achieved by applying a magnetic field perpendicular to the uni-axial
anisotropy field of the free FM [31], as further elaborated in Appendix A4-A.
But it increases the size and power consumption of the device [32]. A recent
work [33] explored the ReLU functionality using Heusler alloy FMs and assuming
only a field-like torque from the SHE layer, which may not give the desired
linear rotation (see appendix: A4-B).

Figure 2: Schematic overview of the simulation setup. (a) Various parameters
of the MTJ and the resistor $R_{1}$ are given as input to the NEGF simulator.
(b) The MTJ resistance is calculated using the NEGF framework, and the voltage
across the MTJ is evaluated from the voltage divider self-consistently. (c) A
schematic of the dependence of the resistance of the MTJ with the FM
polarization angle is shown. (d) The NEGF results are incorporated into the
HSPICE using VerilogA, and the entire circuit is simulated, including the LLGS
equation.

###  II-B Max-Pooling

Max-pooling forms a crucial layer for the CNN, and we discuss in this section
that the proposed ReLU circuit can be appropriately concatenated to perform
the simultaneous local max-pooling function. The max-pooling function
calculates the maximum of the inputs presented [4]. We introduce a competition
among the ReLU circuits to achieve max-pooling with only one winner. We
present a strategy that makes the current input to all the ReLU circuits
(except one) less than zero to have one winner. We will achieve this by
drawing current from the input of the ReLU circuits. The magnitude of the
current drawn depends on the output of other ReLU circuits, thus enabling
competition. At the end of the operation, there will be only one ReLU circuit
with non-zero output. This ReLU output injects a negative current to all other
ReLU devices so that their effective input remains less than zero.
Our design enables the competition through an NMOS transistor connected to a
resistor. The gate terminal of the transistor connects to the output of a
different ReLU circuit. Thus, the current drawn from one ReLU circuit depends
on the output voltage of another ReLU circuit. Figure 1(c) shows the schematic
of the ReLU+max-pooling circuit with two inputs. This circuit design
calculates both the ReLU and max-pooling functions simultaneously. This
strategy makes all outputs zero other than the one corresponding to the
maximum input. We can sum all the outputs and get the correct result without
wondering which device has the maximum input.

Figure 3: Schematic of the device simulation section. The NEGF quantum
transport formalism is self-consistently coupled with the voltage divider
circuit to calculate the MTJ resistance as a function of the FM polarization
angle. Here $H$ is the device Hamiltonian and $G^{n}$ is the electron
correlation matrix, see appendix for more details.

Crossbar arrays [34, 35] can be used to calculate the convolution function,
and the output of the convolution layer is given as input to the proposed
simultaneous ReLU-max-pooling circuit to complete the feature extraction stage
in the CNN as shown in Fig. A1. Non-volatile memory devices such as domain-
wall [36, 3] and skyrmion-based [37, 3] MTJs can be used as synapses in the
crossbar arrays. The input resistance of the proposed network can be tuned (by
changing the dimensions of the SHE layer) as per the loading requirement of
the crossbar array. The proposed ReLU circuit can be directly connected to the
crossbar array without the max-pooling portion to complete the McCulloch-Pitts
neuron[38] that can be used in hidden layers of feed-forward networks [39, 40]
and reservoir computing [41].

TABLE I:  Simulation Parameters. Symbol |  Quantity |  Value
---|---|---
$M_{s}$ |  saturation magnetization |  1150 emu/$cm^{3}$
$H_{k}$ |  anisotropy field |  330 - 3300 Oe[42]
V |  volume of ferromagnet |  1000 $nm^{3}$
$\Delta$ |  thermal stability factor |  4.58 - 45
$\alpha$ |  Gilbert damping |  0.01
L |  length of ferromagnet |  40 nm
W |  width of ferromagnet |  25 nm
$C_{MTJ}$ |  MTJ capacitance |  2.21 fF
$I_{cy}$ |  input current to SHE layer |  35 - 340 $\mu A$
$\theta$ |  spin-hall angle |  0.3 rad [43]
$t_{HM}$ |  thickness of heavy metal |  5 nm
$R_{1}$ |  reference resistor |  698.63 k$\Omega$
$R_{2}$ |  resistance in NMOS current source |  1 - 7 k$\Omega$
$I_{b}$ |  biasing current |  17 - 160 $\mu A$
$V_{DD},V_{SS}$ |  voltage sources |  0.5 V, -0.5 V
$C_{g}$ |  CMOS inverter input capacitance |  0.175 fF
$C_{o}$ |  CMOS inverter output capacitance |  0.305 fF
$\Delta t$ |  simulation time step |  0.5 pS
$\hbar$ |  reduced Plank’s constant |  $1.055\times 10^{-34}$ J s
$k_{B}$ |  Boltzmann constant |  $1.38\times 10^{-16}ergK^{-1}$
T |  temperature |  300 K

##  III Simulation Methods

We show in Fig. 2 the schematic overview of the developed hybrid NEGF-CMOS
simulation framework. The MTJ parameters, such as the effective mass of
electrons in the insulator and the FM, the barrier height of the insulator,
thickness, area, and the resistance value of $R_{1}$ are presented to the NEGF
simulator as shown in Fig. 2(a). The resistance of the MTJ depends on the MTJ
angle and the voltage across MTJ[44, 45], this is shown in Fig. A2a. This
effect is substantial, and using a simple conductance equation to describe the
MTJ ignores this effect. NEGF, on the other hand, includes these effects via
the atomistic level simulation whose results agree with the experimental
data[44, 46, 47]. In circuits employing MTJ, there is usually a change in MTJ
voltage. The combination of these effects stresses the necessity for using
NEGF self-coupled with the circuit simulation to capture the resistance
dependence on the voltage. While simulating CMOS devices such as inverters,
HSPICE is very useful, and integrating NEGF with HSPICE provides a pathway to
design circuits with spintronic and CMOS devices. This integration also allows
access to HSPICE’s analysis tools, such as noise and frequency analysis. This
coupling of NEGF and HSPICE is further helpful in simulating more extensive
networks. The NEGF simulator calculates the resistance of the MTJ for a given
MTJ orientation angle and voltage across it. The voltage across the MTJ
depends on the voltage divider formed by the MTJ and resistor $R_{1}$, so the
NEGF and voltage divider are simulated self consistently as shown in Fig. 3.

This result is coupled to the HSPICE simulator using the VerilogA code. The
HSPICE simulates the entire circuit, including the LLGS equation [48, 49], and
the HSPICE also simulates the approximation of the CMOS inverter pair based on
a 16nm node of the predictive technology model (PTM) [50]. The circuit
parameters used in the simulation are given in Tab.I. We have employed The
Keldysh NEGF formalism[51, 52] coupled with LLGS equations[53, 54, 55] to
describe the magnetization dynamics of the free FM(see appendix:
LABEL:appendix:siml for more details on the simulation methods).

##  IV Results

We show in Fig. 4a the continuous rotation of the magnetization of the free FM
layer with the input current $I_{in}$, the stability of the magnetization is
discussed in appendix A5. The continuous rotation of the magnetization with
$I_{in}$ is achieved when the orthogonal current $I_{cy}$ magnitude is larger
than a critical value (see appendix A4-C). The TMR effect translates the
continuous change in the magnetization to a continuous resistance change of
the MTJ device, as shown in Fig. 4a. The linear region of the MTJ-resistance
around the zero $I_{in}$ has been utilized to obtain the ReLU output as shown
in Fig. 4b by injecting an appropriate bias current $I_{b}$. Various
parameters of the ReLU circuit design are given in Tab. I. The proposed
circuit shown in Fig. 1(b) emulates the ReLU function closely for normalized
inputs of smaller than unity.

(a)

(b)

Figure 4: Magnetization dynamics and the ReLU realization. (a) The
magnetization($m_{y}$) of the free FM layer and the resistance of the MTJ-
device with the current $I_{in}$ along $\hat{x}$-direction, while a constant
critical orthogonal current of 340 $\mu A$ is applied along
$\hat{y}$-direction (see Fig. 1(a)). (b) The normalized output voltage of the
ReLU circuit (see Fig. 1(b)) with normalized input current with $I_{0}=220\mu
A$. The thermal stability factor for the given ReLU is 45.

We evaluate the performance of the ReLU circuit against the thermal stability
factor ($\Delta=\frac{H_{k}M_{s}V}{2k_{B}T}$) of the free FM layer of the MTJ
device. The $\Delta$ factor of the free-FM not only captures the stability of
the magnetization direction against thermal noise but also the extent to which
the spin current changes the magnetization direction. The critical spin
current ($I_{sc}=\frac{4e\alpha k_{B}T\Delta}{\hbar}$) for magnetization
switching is proportional to $\Delta$ [30]. We vary the $\Delta$ of the free
FM layer by changing the anisotropy field, but the same can also be changed by
the MTJ device dimensions. The critical orthogonal current ($I_{cy}$), bias
current ($I_{b}$), and the input current($I_{in}$) required for a continuous
rotation of magnetization decreases with a reduction in $\Delta$. Thus, the
mean static power consumption of the ReLU circuit also decreases with a
reduction in $\Delta$ as shown in Fig. 5a. We also observe that the $\Delta$
reduction does not affect the output settling time ($\approx$4ns) as various
input currents decrease proportionally with $\Delta$.
We show in Fig. 5a the reduction in the output error of the ReLU circuit with
a corresponding increase in $\Delta$. The thermal noise ($\langle
H_{th}^{2}\rangle=\frac{2\alpha k_{B}T}{\gamma M_{s}V}$) stays constant as
$\Delta$ is increased and the weight of the thermal noise decreases with
$H_{eff}$ as $H_{k}$ is increased, since $H_{eff}$ also includes the
contribution from the anisotropy field $H_{k}$. The error in the output is
estimated by performing fifty cycles of Monte Carlo simulations [56] with a
normalized input current ($I_{in}/I_{0}$) of value $0.5$. The decrease in
power consumption and the increase in the error due to the reduction in
$\Delta$ presents an opportunity to optimize the circuit to consume less power
while keeping the error in an acceptable range. It can be inferred from Fig.
5a that the ReLU circuit has an optimal performance of power($0.5-1.2mW$) and
error($<2.2\%$) with $\Delta$ in the range of $15-25$.

Figure 5b shows the boxplot of the error (%) with the $\Delta$. The boxplot
[57] gives profound insights into the error statistics and uses the following
five parameters to summarize the data error. The medians (the red lines)
converging to zero indicate that the average thermal noise is zero. The
reduction in its spread is associated with the variance of the thermal noise
since its effect on free-FM depends on $\Delta$.

(a)

(b)

Figure 5:  Trade-off between power consumption and error percentage. (a)
Static power consumption of the ReLU circuit over the entire input range and
the average of absolute error percentage with normalized input
($I_{in}/I_{0}$) of 0.5. (b) Boxplot of the error percentage of the ReLU
circuit under a thermal noise. The error percentage increases with a decrease
in the thermal stability factor.

We show in Fig. 6a the transient response of the ReLU-max-pooling circuit with
9-inputs corresponding to the typically used $3\times 3$ pooling layer for the
worst-case scenario. In the worst-case scenario, all the inputs are close to
the maximum possible value and the observed settling time of the output is
$9ns$. Initially, all the nine outputs rise to reach the output dictated by
the corresponding inputs. Then quickly, the competition starts, and all the
ReLU circuits inject negative currents into each other through the NMOS
transistors causing the outputs to decrease. As the competition progresses,
one of the outputs (that has the maximum input) becomes the winner and reaches
its steady-state value corresponding to its input while keeping all the other
outputs at zero. It can also be inferred from Fig 6a that the sum of all the
outputs is equal to the maximum output in the steady-state, and it performs
the max-pooling without any knowledge about which of the 9-ReLU devices is
producing it.

(a)

(b)

Figure 6: Max-pooling circuit transient response. (a) The transient response
of the ReLU-max-pooling network with nine inputs. The inputs are chosen to
demonstrate the worst-case scenario. (b) The output of the ReLU-max-pooling
network with nine inputs, the x-axis is the maximum input current out of all
the nine inputs normalized to $I_{0}=220\mu A$.

Figure 6b shows the ReLU-max-pooling network results for $200$ Monte Carlo
simulations. The value of the input currents is varied using the Monte Carlo
method. The inputs are selected such that the maximum input is varied over the
entire input range. The results from Fig. 6b closely resemble the ReLU
function while simultaneously calculating the max pooling.
We show in Fig. 7a the ReLU-max-pooling network’s energy consumption to reach
the steady-state in the worst-case scenario. The worst-case energy consumption
decreases as $\Delta$ decreases. This decrease in energy is due to the
reduction of energy requirement for each ReLU circuit and the reduction in the
competition circuit energy that is caused by a decrease in the NMOS device
dimensions and resistance $R_{2}$ due to the decrease in $I_{in}$ required.

(a)

(b)

(c)

Figure 7: Energy requirements. (a) Worst-case energy requirement for the
output to settle for a 9-input ReLU-max-pooling network. (b) Static power
consumption of a 9-input ReLU-max-pooling circuit averaged over the entire
range of inputs and the average of absolute error percentage of the network
for a normalized input of 0.5. (c) Boxplot of the error percentage of a
9-input network (see Fig. 1(c)) under thermal noise.

Figure 7b shows the static power consumption and percentage error of the
network. The static power consumption decreases, and the percentage error
increases with a reduction in $\Delta$. This behavior shares the same
rationale as that of a single ReLU device discussed earlier. As $\Delta$
increases, there is a linear increase in power consumption and a sharp fall in
error ($\%$). It opens up a possibility for the optimization of the power
consumption ($120\mu W$) with an acceptable percentage error ($2.3\%$) at
$\Delta\approx 15$.
The above discussion suggests that a room temperature unstable free FM
($\Delta<40)$ based MTJ can be utilized to design an energy-efficient ReLU-
max-pooling network. We show in Fig. 7c the boxplot of the error percentage of
the 9-input network. Here the medians converge to zero, indicating that the
average thermal noise is zero. As $\Delta$ increases, the effect of thermal
noise on the free-FM decreases, causing a decline in the spread of the
boxplot.

##  V Conclusion

In this paper, we proposed a circuit design for calculating crucial
neuromorphic functions (ReLU and max-pooling) based on a continuous rotation
of magnetization of the free FM of an MTJ through orthogonal spin current
injection. Using our developed simulation platform based on NEGF and LLGS
equations coupled with the HSPICE circuit simulator, we showed the optimal
range of performance (power and error) against the thermal stability factor of
the free FM of the MTJ device. We demonstrated that the designed circuit is
robust against thermal noise, while consuming $0.56\mu W$ of power for ReLU
functionality and requires $3.4pJ$ of energy for 9-input ReLU-max-pooling
computation in the worst-case scenario, and consumes 120 $\mu W$ of static
power in a typical scenario with a percentage error of 2.3% at $\Delta=15$.
Our work also opens the possibility of using an unstable FM ($\Delta=15$) for
an energy-efficient ReLU-max-pooling circuit with reliable output.

## Acknowledgements

The author BM acknowledges the Visvesvaraya Ph.D Scheme of the Ministry of
Electronics and Information Technology (MEITY), Government of India,
implemented by Digital India Corporation (formerly Media Lab Asia), and also
acknowledges the support by the Science and Engineering Research Board (SERB),
Government of India, Grant No. STR/2019/000030 and Grant No. MTR/2021/000388.
We also acknowledge the Ministry of Human Resource Development (MHRD),
Government of India, Grant No. STARS/APR2019/NS/226/FS under the STARS scheme.
The author AS acknowledges the support of ISIRD phase-1 project of IIT Ropar.

## APPENDIX

##  A1 CNN architecture

An example of CNN architecture using crossbar arrays and ReLU-max pooling
circuits is given in Fig. A1. Here The combination of the crossbar array and
ReLU-max pooling circuit is repeated multiple times (depending on the
complexity of the problem) to complete the feature extraction stage of CNN,
whose result is then given to the classification stage to complete the CNN
network.

##  A2 MTJ resistance dependence on bias voltage

The resistance of the MTJ depends on the MTJ angle and the voltage across MTJ
as shown in Fig. A2a. The voltage change with respect to the input current in
the proposed ReLU circuit is shown in Fig.A2b.

Figure A1: An example of CNN architecture using crossbar arrays and ReLU-max
pooling circuits. The combination of the crossbar array and ReLU-max pooling
circuit is repeated multiple times (depending on the complexity of the
problem) to complete the feature extraction stage of CNN, whose result is then
given to the classification stage to complete the CNN network.

(a)

(b)

Figure A2: (a) Change in resistance of the MTJ due to variation in voltage for
different MTJ angles ($\theta$). (b) The voltage across the MTJ with varying
input spin current.

##  A3 Simulation Method

###  A3-A Quantum Transport

The Keldysh NEGF technique is a versatile formalism used to describe quantum
transport across a variety of set ups featuring charge, spin, heat flow [51,
58, 59, 52] and also can be used to describe transport across hybrid quantum
systems featuring the interplay of transport with the physics of dissipation
and incoherent processes involving phonons, photons and also Cooper pairs [60,
61, 62, 63, 64, 65]. We adapt this formalism [52] to simulate the MTJ
structure with MgO sandwiched between free and fixed FM layers. For a tight-
binding Hamiltonian of the CoFeB, the parameters used are exchange splitting
$\delta=2.15eV$, Fermi energy $E_{f}=2.25eV$, effective mass of MgO barrier is
$m_{OX}=0.18m_{e}$ and that of FM contact is $m_{FM}=0.8m_{e}$ , where $m_{e}$
is the free electron mass. The barrier height of the CoFeB-MgO interface is
$U_{B}=0.76eV$ above the Fermi energy [66, 46], and we have adjusted the width
of the MgO layer to 3nm so that spin-transfer torque has no impact on the
magnetization dynamics of the free FM layer.

We start by describing the Green’s function matrix $[G(E)]$

$[G(E)]=[EI-H-\Sigma]^{-1}$ (1)

$[\Sigma]=[\Sigma_{T}]+[\Sigma_{B}]$ (2)

Where [H] is the device Hamiltonian, $[H]=[H_{0}]+[U]$, comprising device
tight-binding matrix $[H_{0}]$ and the Coulomb charging matrix $[U]$, and [I]
is the identity matrix, $E$ is the energy variable. The self-energy matrix
$[\Sigma]$ is the sum of the top $[\Sigma_{T}]$ and bottom $[\Sigma_{B}]$ FM
layer’s self-energy matrices. The charging matrix $[U]$ is calculated self-
consistently using Poisson’s equation.

The electron correlation matrix $[G^{n}(E)]$ and the in-scattering function
$[\Sigma^{in}(E)]$ are given by

$[G^{n}]=\int dE[G(E)][\Sigma^{in}(E)][G(E)]^{\dagger}$ (3)

$[\Sigma^{in}(E)]=[\Gamma_{T}(E)]f_{T}(E)+[\Gamma_{B}(E)]f_{B}(E)$ (4)

where, $[\Gamma_{T}(E)]=i([\Sigma_{T}(E)]-[\Sigma_{T}(E)]^{\dagger})$ and
$[\Gamma_{B}(E)]=i([\Sigma_{B}(E)]-[\Sigma_{B}(E)]^{\dagger})$ represent spin-
dependent broadening matrices of the top and bottom contacts.

The quantum transport segment [51, 52] culminates with the calculation of
currents via the current operator that describes the current flow between
sites $i$ and $i+1$ is given by

$I_{op}=\frac{i}{\hbar}(H_{i,i+1}G^{n}_{i+1,i}-H_{i+1,i}G^{n}_{i,i+1})$ (5)

This current operator $\hat{I}_{op}$ is $2\times 2$ matrix in the spin space
of a lattice point, using which the charge current can be evaluated as

$I=q\int Real[Trace(\hat{I}_{op})]dE,$ (6)

where $q$ is the quantum of electronic charge. We also consider the effect of
MTJ capacitance since it dominates the CMOS inverter pair capacitance.

###  A3-B Magnetization dynamics

We utilize the LLGS equation [53, 54] to describe the magnetization dynamics
of the free FM. The LLGS equation is given by

$(\frac{1+\alpha^{2}}{\gamma
H_{k}})\frac{d\hat{m}}{dt}=-\hat{m}\times\vec{h}_{eff}-\alpha\hat{m}\times\hat{m}\times\vec{h}_{eff}\\\
-\hat{m}\times\hat{m}\times\vec{i}_{s}+\alpha\hat{m}\times\vec{i}_{s},$ (7)

where $\hat{m}$ is the unit vector along the direction of magnetization of the
free magnet, $\gamma$ is the gyromagnetic ratio, $\alpha$ is the Gilbert
damping parameter, $\vec{h}_{eff}=\frac{\vec{H}_{eff}}{H_{k}}$ is the reduced
effective field and $\vec{i}_{s}=\frac{\hbar\vec{I}_{s}}{2qM_{s}VH_{k}}$ is
the normalized spin current. The term $\vec{H}_{eff}$ includes the
contribution of the anisotropy field ($H_{k}$), the applied magnetic field
($H_{app}$) and the thermal noise ($H_{th}$). The thermal noise is given by
$\langle H_{th}^{2}\rangle=\frac{2\alpha k_{B}T}{\gamma M_{s}V}$ and
$\langle\rangle$ represents the ensemble average [67]. The spin current from
the SHE layer is given by [68, 69]

$I_{s}=\theta\frac{l_{FM}}{t_{HM}}I_{c}\times\sigma$ (8)

where $I_{s}$ is the magnitude of the spin current, $\theta$ is the spin Hall
angle, $l_{FM}$ is the length of the FM, $t_{HM}$ is the thickness of the SHE
layer, $I_{c}$ is the charge current, and $\sigma$ is the polarization of the
spin current.

##  A4 Continuous linear rotation in the magnetization

In this section, we explore various possible routes for the continuous
rotation of magnetization. We also examine the feasibility of each circuit
design.

###  A4-A Static magnetic field driven MTJ

Here, we show that linear resistance can be generated from an MTJ by applying
a static magnetic field orthogonal to the uni-axial anisotropy direction of
the free FM.

(a)

(b)

Figure A3: Device schematic. (a) Schematic of an MTJ device. The PMA magnet is
the free FM layer, and the in-plane magnet is the fixed FM layer. The magnetic
field is applied along the $\hat{y}$-direction to get a linear rotation in the
magnetization of the free layer. (b) Variation of the MTJ resistance with the
normalized magnetic field ($H/H_{k}$) with $V_{MTJ}=0.33V$. Inset shows the
variation of the magnetization ($m_{y}$) of the free FM layer with an applied
magnetic field.

Fig. A3(a) shows the schematic of MTJ with PMA FM as free layer and in-plane
magnet as fixed FM layer. To generate linear rotation in magnetization,
magnetic field is applied in the y-direction. Figure A3(b) shows the
simulation results of the MTJ. LLGS and NEGF formalism are employed in the
simulation of the device setup. This method of using magnetic field is not
recommended because it increases the size and power consumption of the device
[32].

###  A4-B Field-like spin torque driven MTJ

The LLGS equation describes the magnetization dynamics of the free FM.
Equation 9 shows the LLGS equation.

$(\frac{1+\alpha^{2}}{\gamma
H_{k}})\frac{d\hat{m}}{dt}=-\hat{m}\times\vec{h}_{eff}-\alpha\hat{m}\times\hat{m}\times\vec{h}_{eff}\\\
-\hat{m}\times\hat{m}\times\vec{i}_{s}+\alpha\hat{m}\times\vec{i}_{s}$ (9)

(a)

(b)

Figure A4: Magnetization dynamics of the free FM layer for $I_{cy}=340\mu A$
and (a) $I_{in}=266\mu A$ (b) $I_{in}=-266\mu A$ .

The spin current can be resolved as

$\vec{i}_{s}=i_{s,m}\hat{m}+i_{s,\parallel}\hat{N}+i_{s,\perp}\hat{N}\times\hat{m},$
(10)

where $\hat{N}$ is the unit vector along some direction, using (9) and (10)
LLGS equation can be modified as

$(\frac{1+\alpha^{2}}{\gamma
H_{k}})\frac{d\hat{m}}{dt}=-\hat{m}\times\vec{h}_{eff}-\alpha\hat{m}\times\hat{m}\times\vec{h}_{eff}\\\
-\hat{m}\times\hat{m}\times(i_{s,\parallel}+\alpha i_{s,\perp})\hat{N}\\\
-\hat{m}\times(i_{s,\perp}-\alpha i_{s,\parallel})\hat{N}\\\ $ (11)

From (11) it can be inferred that the spin current $(i_{s,\parallel}+\alpha
i_{s,\perp})$ (the Slonczewski term) acts as a damping/anti-damping term and
$(i_{s,\perp}-\alpha i_{s,\parallel})$ acts as a precession term with field-
like behaviour along $\hat{N}$. To generate linear rotation in magnetization
of the free FM and linear resistance change in the MTJ , the magnitude of the
field-like term needs to be increased with $\hat{N}$ along $\hat{y}$. A fixed
$i_{s,\perp}$ cannot be applied to the FM using the spin hall effect as its
magnitude ($\hat{N}\times\hat{m}$) depends on the magnetization of the FM
($\hat{m}$) and changes continuously as the magnetization is rotated. The
field-like term can also be changed via $i_{s,\parallel}$ whose direction
(along $\hat{N}$) is fixed, but it requires 1/$\alpha$ ($\approx 100$) times
more current. And $i_{s,\parallel}$ also affects the magnetization through the
damping term. So, this method is not feasible for generating continuous
rotation of magnetization of free-FM.

###  A4-C Orthogonal spin current driven MTJ

Our design uses orthogonal spin currents to generate the continuous rotation
of magnetization in the free FM layer. These spin currents are produced by
injecting two currents into the heavy metal layer.

Figure A4 shows the magnetization dynamics of the rotation. The dynamics are
examined for $I_{in}=266\mu A$ and $-266\mu A$ while the orthogonal
current($I_{cy}$) is kept at $340\mu A$. The settling time measured is
$0.2ns$.

Figure A5: Magnetization($m_{y}$) of free FM with varying $I_{in}$, for
multiple $I_{cy}$.

One of the orthogonal currents is fixed ($I_{cy}$) to get the continuous
rotation in magnetization, while the other current($I_{in}$) is varied. The
$I_{cy}$ current needs to cross a threshold to generate the continuous
rotation. Fig. A5 shows that $I_{cy}$ needs to be greater than $334\mu A$ to
achieve the continuous rotation in magnetization of the free FM layer for
$\Delta=45$.

##  A5 Stability of Magnetization

In this section, we explore different ways to validate the stability of
magnetization.

###  A5-A Energy profile

The stability of a system can be usually found from the energy profile. But,
in the presence of the spin current, the energy profile
($E=-\vec{M}\cdot\vec{H}_{eff}$) of the free ferromagnet (FM) does not change
significantly due to the non-conservative nature of the Slonzewski spin torque
[70].

In the presence of spin current, the LLGS equation is given by

$(\frac{1+\alpha^{2}}{\gamma})\frac{d\hat{m}}{dt}=-\hat{m}\times(\vec{H}_{eff}-\alpha\beta\vec{I}_{s})-\hat{m}\times\hat{m}\times(\alpha\vec{H}_{eff}+\beta\vec{I}_{s})$
(12)

Here, $\beta=\frac{\hbar}{2qM_{s}V}$, the spin current $I_{s}$, and $\alpha$
is the damping factor. For small value of $\alpha$,
$\alpha\beta\vec{I}_{s}<<\vec{H}_{eff}$ the LLGS equation can be rewritten as

$(\frac{1+\alpha^{2}}{\gamma})\frac{d\hat{m}}{dt}=-\hat{m}\times\vec{H}_{eff}-\hat{m}\times\hat{m}\times(\alpha\vec{H}_{eff}+\beta\vec{I}_{s})$
(13)

It can be seen from Equation 2 that the spin current in the LLGS equation act
either as a damping torque or an anti-damping torque. Being a non-conservative
term [70], the spin torque can not affect the energy profile of the FM. One
may draw an analogy with the friction of motion, which may either opposes the
motion (damping term: always the case for friction) or may hypothetically
support the motion (anti-damping), but it can not change the conservation
energy profile (e.g., gravitational potential energy) of the system. Hence, we
can not employ the energy profile of the ferromagnet in the presence of
orthogonal spin currents (in our case) to analyze the stability of the
magnetization.

(a)

(b)

(c)

(d)

(e)

(f)

Figure A6: (a), (b) Magnetization dynamics of free-FM of the MTJ with an
applied spin current of 510$\mu$A with x-polarization and 400$\mu$A with
y-polarization. It can be seen here that the magnetization is stable after it
reaches its final value of $m_{x}=0.7869,m_{y}=0.6171,m_{z}=0.00$ with
different initial magnetizations of -z and +z. (c), (d), (e) and (f) Noise
analysis of magnetization dynamics of the free FM with 1000 simulations of
LLGS equation under the influence of thermal noise, having a thermal stability
factor($\Delta$) of $45$. (c) Variation in $m_{x}$ magnetization due to
thermal noise. (d) Variation in $m_{y}$ magnetization due to thermal noise.
(e) Variation in $m_{z}$ magnetization due to thermal noise. (f) Boxplot of
error $\%$ of $m_{x}$ and $m_{y}$ magnetization.

###  A5-B LLG simulation

The stability of the magnetization in the presence of the spin current can be
analyzed via evaluation of the magnetization from the different initial
configurations and/or stochastic evaluation of the magnetization with thermal
noise ($\langle H_{th}^{2}\rangle=\frac{2\alpha k_{B}T}{\gamma M_{s}V}$). We
show in Fig. A6a&(b) the magnetization dynamics of the FM in the presence of
orthogonal spin currents 510$\mu$A, 400$\mu$A with x-and y-polarizations
respectively, producing a stable magnetization orientation of
$m_{x}=0.7869,m_{y}=0.6171,m_{z}=0.00$. It can be observed (Fig. A6a&(b)) the
different initial magnetization of the FM stabilizes to the same final value
for given spin current inputs. Figures A6c,A6d and A6e show the effect of
thermal noise on $m_{x},m_{y}$, and $m_{z}$ magnetizations respectively. The
error percentage of magnetization is shown in Fig. A6f, it can be inferred
that the thermal noise has minimal effect on the magnetization with error
$\leq 1\%$, leading to a stable magnetization rotation in the presence of the
spin currents. The error is estimated by performing 1000 stochastic
simulations.

###  A5-C Net torque on the FM

The stability of magnetization can be analyzed from the plots of torque vs
polar angle ($\theta$) and azimuthal angle ($\phi$) of the magnetization.

(a)

(b)

(c)

Figure A7: (a) Norm of the net torque acting on the free-FM in the presence of
the spin current of $500\mu A$ with z-polarization. The net torque vanishes at
$\theta=0$ and at $\theta=\pi$. (b) The norm of the net torque with $\theta$
at $\phi=0$ is plotted to clearly spot the vanishing torque points. (c) The
in-plane torque with $\theta$ for the vanishing torque points at $\phi=0$.

(a)

(b)

(c)

Figure A8: (a) Norm of the net torque acting on the free-FM in the presence of
the spin current of $510\mu A$ and $600\mu A$ with x & y-polarization,
respectively. The net torque vanishes at $m_{1}$ and $m_{2}$ points. (b) The
norm of the net torque with $\phi$ at $\theta=\pi/2$ is plotted to clearly
spot the vanishing torque points. (c) The in-plane torque with $\phi$ for the
vanishing torque points at $\theta=\pi/2$.

The net torque on the free FM is given by:

$\vec{\tau}=-\hat{m}\times\vec{H}_{eff}-\alpha\hat{m}\times\hat{m}\times\vec{H}_{eff}-\beta\hat{m}\times\hat{m}\times\vec{I}_{s}+\alpha\beta\hat{m}\times\vec{I}_{s}$
(14)

The net torque($\tau$) can be decomposed into in-plane torque($\tau_{ip}$) and
out-of-plane torque ($\tau_{op}$). The $\tau_{ip}$ can be damping (positive)
or anti-damping (negative). The $\tau_{ip}$ is given by

$\tau_{ip}=\hat{i}_{p}\cdot\vec{\tau}$ (15)

where, $\hat{i}_{p}$ is given by

$\hat{i}_{p}=\frac{-\hat{m}\times\hat{m}\times\hat{m}_{v}}{|-\hat{m}\times\hat{m}\times\hat{m}_{v}|}\\\
$ (16)

Here, $\hat{m}_{v}$ is the vanishing torque point($\vec{\tau}=0$)[71], which
can be extracted from plots of torque vs $\theta$ and $\phi$ of the
magnetization. If the $\tau_{ip}$ is damping(positive), that means the
$\hat{m}_{v}$ is a stable point, and if it is anti-damping(negative), then
$\hat{m}_{v}$ is not a stable point in the presence of the applied spin
current. We first demonstrate the effectiveness of this method on regular
switching(anti-parallel to parallel, i.e., PMA injected with z-polarized spin
current or vice versa) and adapt the same to explain the magnetization
rotation via orthogonal spin current injection considered in the paper.

### Regular switching

Consider a PMA ferromagnet injected with a spin current of $500\mu A$ with
z-polarization, higher than the critical spin current for switching. It is
well known that the resulting magnetization will be aligned in the
+z-direction. So, in the presence of spin current, +z is the stable point,
whereas -z is not a stable point.

Figure A7a shows the norm of the total torque acting on the free FM. It shows
two regions ($\theta=0,\pi$ pointing to +z and -z directions, respectively)
where the torque becomes zero which can be clearly seen from Fig. A7b. So,
even in the presence of spin current, the net torque vanishes at the -z
direction, this means that the net torque alone cannot fully explain the
stability of the magnetization.

To fully verify the stability of the magnetization point, we can look at the
in-plane torque($\tau_{ip}$). Figure A7c shows the in-plane torque for both of
the vanishing torque points. It can see that for $\hat{m}_{v}=$[0 0 1], the
$\tau_{ip}$ is always positive(damping), implying a stable magnetization
point, whereas for $\hat{m}_{v}=$[0 0 -1], the $\tau_{ip}$ is always
negative(anti-damping) implying not a stable magnetization point in the
presence of the spin current.

### Orthogonal spin current injected free FM

The stability of the magnetization considered in this paper can be explored in
a similar fashion. The spin currents of $510\mu A$ and $600\mu A$ with x &
y-polarization, respectively, are injected into the PMA free-FM, and the
expected stable point for this condition is $m=$[0.6477 0.7619 0]. We show in
Fig. A8a the net torque on the free-FM. It can be seen that net torque
vanishes at two points($m_{1}=$[0.6476 0.7619 0] and $m_{2}=$[-0.6476 -0.7619
0]), and both the torque vanishing points appear at $\theta=\pi/2$, these
points can be clearly seen from Fig. A8b.

We show the in-plane torque in Fig. A8c to inspect the stability of these
torque vanishing points. The $\tau_{ip}$ is always positive for $m_{1}$
implying the stable point, and always negative for $m_{2}$, implying not a
stable point in the presence of the orthogonal spin currents. We believe that
the thermal stability of the magnetization may be correlated to the area under
the in-plane torque curve, which will be addressed in our future works.

## References

  * [1] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, “Physics for neuromorphic computing,” _Nature Reviews Physics_ , vol. 2, no. 9, pp. 499–510, 2020.
  * [2] C. D. Schuman, T. E. Potok, R. M. Patton, J. D. Birdwell, M. E. Dean, G. S. Rose, and J. S. Plank, “A survey of neuromorphic computing and neural networks in hardware,” _arXiv preprint arXiv:1705.06963_ , 2017.
  * [3] J. Grollier, D. Querlioz, K. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, “Neuromorphic spintronics,” _Nature electronics_ , vol. 3, no. 7, pp. 360–370, 2020.
  * [4] K. O’Shea and R. Nash, “An introduction to convolutional neural networks,” _arXiv preprint arXiv:1511.08458_ , 2015.
  * [5] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [6] S. Khan, H. Rahmani, S. A. A. Shah, and M. Bennamoun, “A guide to convolutional neural networks for computer vision,” _Synthesis Lectures on Computer Vision_ , vol. 8, no. 1, pp. 1–207, 2018.
  * [7] L. Xie and A. Yuille, “Genetic cnn,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 1379–1388.
  * [8] S. Albawi, T. A. Mohammed, and S. Al-Zawi, “Understanding of a convolutional neural network,” in _2017 International Conference on Engineering and Technology (ICET)_. Ieee, 2017, pp. 1–6.
  * [9] S. Sharma, S. Sharma, and A. Athaiya, “Activation functions in neural networks,” _towards data science_ , vol. 6, no. 12, pp. 310–316, 2017.
  * [10] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT press, 2016.
  * [11] Y.-L. Boureau, J. Ponce, and Y. LeCun, “A theoretical analysis of feature pooling in visual recognition,” in _Proceedings of the 27th international conference on machine learning (ICML-10)_ , 2010, pp. 111–118.
  * [12] X. Glorot, A. Bordes, and Y. Bengio, “Deep sparse rectifier neural networks,” in _Proceedings of the fourteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 2011, pp. 315–323.
  * [13] C.-H. Chang, H.-Y. Kao, and S.-H. Huang, “Hardware implementation for multiple activation functions,” in _2019 IEEE International Conference on Consumer Electronics-Taiwan (ICCE-TW)_. IEEE, 2019, pp. 1–2.
  * [14] C. Geng, Q. Sun, and S. Nakatake, “An analog cmos implementation for multi-layer perceptron with relu activation,” in _2020 9th International conference on modern circuits and systems technologies (MOCAST)_. IEEE, 2020, pp. 1–6.
  * [15] P. Priyanka, G. Nisarga, and S. Raghuram, “Cmos implementations of rectified linear activation function,” in _International Symposium on VLSI Design and Test_. Springer, 2018, pp. 121–129.
  * [16] S. Furber, “Large-scale neuromorphic computing systems,” _Journal of neural engineering_ , vol. 13, no. 5, p. 051001, 2016.
  * [17] G. Indiveri, B. Linares-Barranco, T. J. Hamilton, A. v. Schaik, R. Etienne-Cummings, T. Delbruck, S.-C. Liu, P. Dudek, P. Häfliger, S. Renaud _et al._ , “Neuromorphic silicon neuron circuits,” _Frontiers in neuroscience_ , vol. 5, p. 73, 2011.
  * [18] A. Parihar, N. Shukla, M. Jerry, S. Datta, and A. Raychowdhury, “Computing with dynamical systems based on insulator-metal-transition oscillators,” _Nanophotonics_ , vol. 6, no. 3, pp. 601–611, 2017.
  * [19] J. Cai, B. Fang, L. Zhang, W. Lv, B. Zhang, T. Zhou, G. Finocchio, and Z. Zeng, “Voltage-controlled spintronic stochastic neuron based on a magnetic tunnel junction,” _Physical Review Applied_ , vol. 11, no. 3, p. 034015, 2019.
  * [20] A. Sengupta and K. Roy, “Spin-transfer torque magnetic neuron for low power neuromorphic computing,” in _2015 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2015, pp. 1–7.
  * [21] J. Cai, L. Zhang, B. Fang, W. Lv, B. Zhang, G. Finocchio, R. Xiong, S. Liang, and Z. Zeng, “Sparse neuromorphic computing based on spin-torque diodes,” _Applied Physics Letters_ , vol. 114, no. 19, p. 192402, 2019.
  * [22] E. Raimondo, A. Giordano, A. Grimaldi, V. Puliafito, M. Carpentieri, Z. Zeng, R. Tomasello, and G. Finocchio, “Reliability of neural networks based on spintronic neurons,” _IEEE Magnetics Letters_ , vol. 12, pp. 1–5, 2021.
  * [23] V. B. Desai, D. Kaushik, J. Sharda, and D. Bhowmik, “On-chip learning of a domain-wall-synapse-crossbar-array-based convolutional neural network,” _Neuromorphic Computing and Engineering_ , vol. 2, no. 2, p. 024006, 2022\.
  * [24] W. A. Borders, H. Akima, S. Fukami, S. Moriya, S. Kurihara, Y. Horio, S. Sato, and H. Ohno, “Analogue spin–orbit torque device for artificial-neural-network-based associative memory operation,” _Applied physics express_ , vol. 10, no. 1, p. 013007, 2016.
  * [25] U. Sahu, N. Sisodia, P. K. Muduli, and D. Bhowmik, “Ferrimagnetic synapse devices for fast and energy-efficient on-chip learning on crossbar-array-based neural networks (a device-circuit-system costudy),” _IEEE Transactions on Electron Devices_ , vol. 69, no. 4, pp. 1713–1720, 2022\.
  * [26] U. Singh, N. Garg, S. Kumar, P. K. Muduli, and D. Bhowmik, “Learning of classification tasks with an array of uniform-mode spin hall nano-oscillators,” _AIP Advances_ , vol. 11, no. 4, p. 045117, 2021.
  * [27] H. Abbas, Y. Abbas, S. N. Truong, K.-S. Min, M. R. Park, J. Cho, T.-S. Yoon, and C. J. Kang, “A memristor crossbar array of titanium oxide for non-volatile memory and neuromorphic applications,” _Semiconductor Science and Technology_ , vol. 32, no. 6, p. 065014, 2017.
  * [28] M. Hu, J. P. Strachan, Z. Li, E. M. Grafals, N. Davila, C. Graves, S. Lam, N. Ge, J. J. Yang, and R. S. Williams, “Dot-product engine for neuromorphic computing: Programming 1t1m crossbar to accelerate matrix-vector multiplication,” in _2016 53nd acm/edac/ieee design automation conference (dac)_. IEEE, 2016, pp. 1–6.
  * [29] Y. Li and K.-W. Ang, “Hardware implementation of neuromorphic computing using large-scale memristor crossbar arrays,” _Advanced Intelligent Systems_ , vol. 3, no. 1, p. 2000137, 2021.
  * [30] A. Sharma, A. A. Tulapurkar, and B. Muralidharan, “Proposal for energy efficient spin transfer torque-magnetoresistive random access memory device,” _Journal of Applied Physics_ , vol. 129, no. 23, p. 233901, 2021\.
  * [31] A. Sharma, A. Tulapurkar, and B. Muralidharan, “Ultrasensitive nanoscale magnetic-field sensors based on resonant spin filtering,” _IEEE Transactions on Electron Devices_ , vol. 63, no. 11, pp. 4527–4534, 2016.
  * [32] D. Apalkov, B. Dieny, and J. Slaughter, “Magnetoresistive random access memory,” _Proceedings of the IEEE_ , vol. 104, no. 10, pp. 1796–1830, 2016\.
  * [33] A. W. Stephan and S. J. Koester, “Spin hall mtj devices for advanced neuromorphic functions,” _IEEE Transactions on Electron Devices_ , vol. 67, no. 2, pp. 487–492, 2020.
  * [34] S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, and W. Lu, “Nanoscale memristor device as synapse in neuromorphic systems,” _Nano letters_ , vol. 10, no. 4, pp. 1297–1301, 2010.
  * [35] D. Fan, M. Sharad, and K. Roy, “Design and synthesis of ultralow energy spin-memristor threshold logic,” _IEEE Transactions on Nanotechnology_ , vol. 13, no. 3, pp. 574–583, 2014.
  * [36] S. A. Siddiqui, S. Dutta, A. Tang, L. Liu, C. A. Ross, and M. A. Baldo, “Magnetic domain wall based synaptic and activation function generator for neuromorphic accelerators,” _Nano letters_ , vol. 20, no. 2, pp. 1033–1040, 2019.
  * [37] Y. Huang, W. Kang, X. Zhang, Y. Zhou, and W. Zhao, “Magnetic skyrmion-based synaptic devices,” _Nanotechnology_ , vol. 28, no. 8, p. 08LT02, 2017.
  * [38] B. Goertzel, _The structure of intelligence: A new mathematical model of mind_. Springer Science & Business Media, 2013.
  * [39] I. W. Sandberg, J. T. Lo, C. L. Fancourt, J. C. Principe, S. Katagiri, and S. Haykin, _Nonlinear dynamical systems: feedforward neural network perspectives_. John Wiley & Sons, 2001, vol. 21.
  * [40] D. Svozil, V. Kvasnicka, and J. Pospichal, “Introduction to multi-layer feed-forward neural networks,” _Chemometrics and intelligent laboratory systems_ , vol. 39, no. 1, pp. 43–62, 1997.
  * [41] K. Nakajima and I. Fischer, _Reservoir Computing_. Springer, 2021.
  * [42] M. Gajek, J. Nowak, J. Sun, P. Trouilloud, E. O’sullivan, D. Abraham, M. Gaidis, G. Hu, S. Brown, Y. Zhu _et al._ , “Spin torque switching of 20 nm magnetic tunnel junctions with perpendicular anisotropy,” _Applied Physics Letters_ , vol. 100, no. 13, p. 132408, 2012.
  * [43] C.-F. Pai, L. Liu, Y. Li, H. Tseng, D. Ralph, and R. Buhrman, “Spin transfer torque devices utilizing the giant spin hall effect of tungsten,” _Applied Physics Letters_ , vol. 101, no. 12, p. 122404, 2012.
  * [44] J. C. Sankey, Y.-T. Cui, J. Z. Sun, J. C. Slonczewski, R. A. Buhrman, and D. C. Ralph, “Measurement of the spin-transfer-torque vector in magnetic tunnel junctions,” _Nature Physics_ , vol. 4, no. 1, pp. 67–71, 2008.
  * [45] C. Wang, Y.-T. Cui, J. Sun, J. Katine, R. Buhrman, and D. Ralph, “Bias and angular dependence of spin-transfer torque in magnetic tunnel junctions,” _Physical Review B_ , vol. 79, no. 22, p. 224416, 2009.
  * [46] H. Kubota, A. Fukushima, K. Yakushiji, T. Nagahama, S. Yuasa, K. Ando, H. Maehara, Y. Nagamine, K. Tsunekawa, D. D. Djayaprawira _et al._ , “Quantitative measurement of voltage dependence of spin-transfer torque in mgo-based magnetic tunnel junctions,” _Nature Physics_ , vol. 4, no. 1, pp. 37–41, 2008.
  * [47] G.-X. Miao, K. B. Chetry, A. Gupta, W. H. Butler, K. Tsunekawa, D. Djayaprawira, and G. Xiao, “Inelastic tunneling spectroscopy of magnetic tunnel junctions based on cofeb/ mgo/ cofeb with mg insertion layer,” _Journal of applied physics_ , vol. 99, no. 8, p. 08T305, 2006.
  * [48] G. D. Panagopoulos, C. Augustine, and K. Roy, “Physics-based spice-compatible compact model for simulating hybrid mtj/cmos circuits,” _IEEE Transactions on Electron Devices_ , vol. 60, no. 9, pp. 2808–2814, 2013.
  * [49] J. Z. Sun, “Spin-current interaction with a monodomain magnetic body: A model study,” _Physical Review B_ , vol. 62, no. 1, p. 570, 2000.
  * [50] “Predictive technology model (ptm),” http://ptm.asu.edu/, (Accessed on 11/23/2021).
  * [51] S. Datta, _Electronic transport in mesoscopic systems_. Cambridge university press, 1997.
  * [52] A. Sharma, A. A. Tulapurkar, and B. Muralidharan, “Resonant spin-transfer-torque nano-oscillators,” _Physical Review Applied_ , vol. 8, no. 6, p. 064014, 2017.
  * [53] J. C. Slonczewski, “Current-driven excitation of magnetic multilayers,” _Journal of Magnetism and Magnetic Materials_ , vol. 159, no. 1-2, pp. L1–L7, 1996.
  * [54] A. Brataas, A. D. Kent, and H. Ohno, “Current-induced torques in magnetic materials,” _Nature materials_ , vol. 11, no. 5, pp. 372–381, 2012.
  * [55] A. Sharma, A. A. Tulapurkar, and B. Muralidharan, “Role of phase breaking processes on resonant spin transfer torque nano-oscillators,” _AIP Advances_ , vol. 8, no. 5, p. 055913, 2018.
  * [56] J. Hammersley, _Monte carlo methods_. Springer Science & Business Media, 2013.
  * [57] M. Hubert and E. Vandervieren, “An adjusted boxplot for skewed distributions,” _Computational statistics & data analysis_, vol. 52, no. 12, pp. 5186–5201, 2008.
  * [58] D. Datta, “Modeling of spin transport in mtj devices,” Ph.D. dissertation, Purdue University, 2012.
  * [59] S. Datta, _Lessons from Nanoelectronics: A New Perspective on Transport—Part B: Quantum Transport_. World Scientific, 2018.
  * [60] A. Singha and B. Muralidharan, “Performance analysis of nanostructured peltier coolers,” _Journal of Applied Physics_ , vol. 124, no. 14, p. 144901, 2018\. [Online]. Available: https://doi.org/10.1063/1.5044254
  * [61] U. Aeberhard, “Quantum-kinetic theory of photocurrent generation via direct and phonon-mediated optical transitions,” _Phys. Rev. B_ , vol. 84, p. 035454, Jul 2011. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.84.035454
  * [62] A. Martín-Rodero, F. J. García-Vidal, and A. Levy Yeyati, “Microscopic theory of josephson mesoscopic constrictions,” _Phys. Rev. Lett._ , vol. 72, pp. 554–557, Jan 1994. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.72.554
  * [63] P. San-Jose, J. Cayao, E. Prada, and R. Aguado, “Multiple andreev reflection and critical current in topological superconducting nanowire junctions,” _New Journal of Physics_ , vol. 15, no. 7, p. 075019, jul 2013. [Online]. Available: https://doi.org/10.1088/1367-2630/15/7/075019
  * [64] P. Sriram, S. S. Kalantre, K. Gharavi, J. Baugh, and B. Muralidharan, “Supercurrent interference in semiconductor nanowire josephson junctions,” _Phys. Rev. B_ , vol. 100, p. 155431, Oct 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.100.155431
  * [65] C. Duse, P. Sriram, K. Gharavi, J. Baugh, and B. Muralidharan, “Role of dephasing on the conductance signatures of majorana zero modes,” _Journal of Physics: Condensed Matter_ , vol. 33, no. 36, p. 365301, jul 2021\. [Online]. Available: https://doi.org/10.1088/1361-648x/ac0d16
  * [66] D. Datta, B. Behin-Aein, S. Datta, and S. Salahuddin, “Voltage asymmetry of spin-transfer torques,” _IEEE Transactions on Nanotechnology_ , vol. 11, no. 2, pp. 261–272, 2011.
  * [67] J. Z. Sun, T. Kuan, J. Katine, and R. H. Koch, “Spin angular momentum transfer in a current-perpendicular spin-valve nanomagnet,” in _Quantum Sensing and Nanophotonic Devices_ , vol. 5359. International Society for Optics and Photonics, 2004, pp. 445–455.
  * [68] J. Hirsch, “Spin hall effect,” _Physical review letters_ , vol. 83, no. 9, p. 1834, 1999.
  * [69] S. Takahashi and S. Maekawa, “Spin current, spin accumulation and spin hall effect,” _Science and Technology of Advanced Materials_ , vol. 9, no. 1, p. 014105, 2008.
  * [70] S. Salahuddin, D. Datta, and S. Datta, “Spin transfer torque as a non-conservative pseudo-field,” _arXiv preprint arXiv:0811.3472_ , 2008.
  * [71] K. Roy, S. Bandyopadhyay, and J. Atulasimha, “Metastable state in a shape-anisotropic single-domain nanomagnet subjected to spin-transfer-torque,” _Applied Physics Letters_ , vol. 101, no. 16, p. 162405, 2012.