# A High Throughput Generative Vector Autoregression Model for Stochastic
Synapses

T. Hennen 1, A. Elias 2, J. F. Nodin 3, G. Molas 3,4, R. Waser 1, D. J.
Wouters 1, and D. Bedau 2,∗

###### Abstract

## 1

By imitating the synaptic connectivity and plasticity of the brain, emerging
electronic nanodevices offer new opportunities as the building blocks of
neuromorphic systems. One challenge for large-scale simulations of
computational architectures based on emerging devices is to accurately capture
device response, hysteresis, noise, and the covariance structure in the
temporal domain as well as between the different device parameters. We address
this challenge with a high throughput generative model for synaptic arrays
that is based on a recently available type of electrical measurement data for
resistive memory cells. We map this real world data onto a vector
autoregressive stochastic process to accurately reproduce the device
parameters and their cross-correlation structure. While closely matching the
measured data, our model is still very fast; we provide parallelized
implementations for both CPUs and GPUs and demonstrate array sizes above one
billion cells and throughputs exceeding one hundred million weight updates per
second, above the pixel rate of a 30 frames/s 4K video stream.

##  2 Keywords:

Neuromorphic Computing, Machine Learning, Emerging Technologies, Stochastic
Model, Synapse, ReRAM, Julia, GPU

##  3 Introduction

Recent trends in computing hardware have placed increasing emphasis on
neuromorphic architectures implementing machine learning (ML) algorithms
directly in hardware. Such bio-inspired approaches, through in-memory
computation and massive parallelism, excel in new classes of computational
problems and offer promising advantages with respect to power consumption
error resiliency. While CMOS-based neuromorphic computing (NC) implementations
have made substantial progress recently, new materials and physical mechanisms
may ultimately provide better opportunities for energy efficiency and scaling
[1, 2, 3].

Figure 1:  In analogy to biological synapses, two terminal solid state
nanodevices such as ReRAM can store synaptic weights as electrical resistance
states. The devices, consisting simply of patterned metal-insulator-metal
material stacks, have an adjustable resistance level determined by the ionic
configuration inside the insulating layer. This nano-ionic mechanism also
exhibits non-ideal properties such as stochasticity and noise.

A specific functionality required in NC applications is the ability to mimic
synaptic connections and plasticity by allowing the storage of large numbers
of interconnected and continuously adaptable resistance values. Several
candidate memory technologies such as MRAM, ReRAM, PCM, CeRAM, are emerging to
cover this behavior using different physical mechanisms [4, 5, 6, 7]. Among
these, ReRAM is attractive for its simplicity of materials and device
structure, providing the necessary CMOS compatibility and scalability [8].
ReRAM is essentially a two terminal nanoscale electrochemical cell, whose
variable resistance state is based on manipulation of the point defect
configuration in the oxide material (depicted in Fig. 1). This redox-based
switching mechanism is intrinsically analog, allowing stable resistance levels
to be stored and adjusted through application of bipolar voltage stimuli.
However, non-idealities such as stochasticity, nonlinearity, and noise are
prominent features of these devices that critically impact the performance of
neuromorphic systems composed of them [9].

Modern ML models have reached an astonishingly large and ever-increasing size,
with recent examples exceeding a hundred billion weights [10]. Before
comparable neuromorphic hardware using artificial solid-state synapses can
become a reality, large-scale network designs need to first be implemented and
evaluated in computer simulations. Training, validation, and optimization of
such networks is a process that involves a huge number of simulated devices,
voltage pulses, and current readouts. Within this process, it is important to
accurately consider the constraints of the underlying hardware in detail.
Therefore, lightweight, fast, and accurate stochastic simulations of the
individual synaptic devices is a key requirement.

Traditionally, device modelling begins with a physical description of the
materials and processes involved. In the case of ReRAM, the physical situation
is immensely complicated with many degrees of freedom, and accurate modelling
is a wide-scale and ongoing research undertaking. Efforts in this direction
are motivated by advancing an understanding of physical and chemical
dependencies that can in principle inform design choices on physically
justified grounds. In the past decade, many different computational techniques
have been employed to furnish device models, from ab initio density-functional
theory (DFT), molecular dynamics (MD), kinetic Monte Carlo (KMC), finite
element method (FEM), as well as ordinary differential equation (ODE) and
differential algebraic equation (DAE) solvers [11, 12, 13, 14, 15]. The
resulting models exist on a spectrum of physical abstraction, such that the
cost of increasing computational speed is generally a trade-off in physical
accuracy/detail [16].

Device models that naturally encompass stochasticity do so at the cost of
complexity needed to compute the physical scenario in high detail. For
example, atomistic KMC simulates switching processes with atomic precision and
is inherently stochastic but requires hours of computation per cycle even for
small individual cell volumes (e.g. 125 nm2 [17]). At the other end of the
spectrum, dynamic models based on numerical solutions of ODE systems are
designed to run significantly faster while sometimes aiming to remain
physically realistic. However, their higher speed invariably comes at the cost
of approximations, simplifications, and omissions of physical reality.
Typically, device operation is distilled to a dynamical description of one or
two state variables, such as a conducting filament length, radius, or a defect
concentration.

Due in part to ambiguity in their high dimensional parameter space, a given
ODE model encompasses a diverse range of possible cell behaviors and has the
flexibility to approximately match measurement data [18, 19]. However, fitting
the model to data is commonly an ad-hoc, manual, and/or unspecified procedure.
Having dispensed with the atomistic sources of variability, ODE models are
fully deterministic by default. Where stochasticity is required, it is
accounted for by injecting noise into the state variables or parameters of the
model [20, 21, 22]. Due to the unique experimental challenges posed by
electrical measurement of ReRAM, the data used for fitting is not necessarily
statistically sufficient nor measured under relevant electrical conditions and
timescales. While models can be tuned by hand to roughly match the dispersion
observed in a measurement [23, 24], they generally fail to accurately
reproduce the complex statistical properties of actual devices.

The main purpose of ODE device models is to be computationally efficient
enough to support circuit simulation. Still, nonlinear ODE solvers require
many finely spaced timesteps and a considerable amount of total time to
compute dynamical trajectories. Although they have been successfully used to
demonstrate small scale circuitry such as logic elements and small crossbar
arrays [25, 26, 27, 28], benchmarks or indications of run time for ODE-based
simulations have so far not been supplied. With the exception of extremely
small ML model sizes on the order of $10^{3}$ weights or below, demonstrations
of network performance are expected to remain computationally intractable via
conventional circuit simulation.

In this article, we address these device modelling challenges with a new type
of generative model for arrays of artificial synapses. The main objective of
the model is to accurately reproduce the statistical properties of fabricated
devices while remaining computationally lightweight. Starting with newly
available electrical measurement data as an input, this phenomenological model
is systematically fit using a well defined statistical regression analysis.
The exclusive use of easily computable analytical expressions provides close
quantitative agreement with relevant experimental observation. Taking
advantage of parallel resources on a modern CPU and GPU, we demonstrate the
ability to simulate hundreds of millions of synaptic connections with over
$10^{8}$ weight updates per second. With its high throughput and low memory
footprint, the model can be usefully employed to simulate large arrays of
solid-state synapses for investigation of emerging NC concepts on a relevant
scale.

##  4 Method

The basic requirement for an electronic device serving as an artificial
synapse is to moderate the flow of electrical signals through connections in a
network. Left undisturbed, the device ideally maintains a fixed weight, or
dependence between the voltage across the two device terminals, $U$, and the
resulting current through the device, $I$. Further, for learning there must be
some means of affecting the weight in a durable way. ReRAMs are bipolar
devices that have an adjustable (potentially nonlinear) non-volatile
resistance state, which is based on the size and shape of a conducting
filament that partially or fully bridges the insulating gap of the oxide
material. Simplistically, when $U$ exceeds certain threshold levels, the
resistance state begins to transition toward lower or higher values depending
on the voltage polarity, which corresponds to growth and shrinkage of the
conducting filament. When the filament only partially bridges the insulating
gap, conduction may be limited for example by tunneling through a Schottky
barrier of a material interface, leading to a relatively high resistance
levels [29, 8]. As the filament grows and gradually bridges the gap, the
resistance decreases as conduction transitions into the ohmic type.

Figure 2:  Resistance states reached in a synaptic ReRAM device through
application of voltage pulses exhibit a probabilistic dependence on past
states, leading to long range correlations that also involve other parameters
such as the voltage thresholds required for switching. Starting with
effectively infinite state possibilities, represented by the three cells on
the left, an applied voltage pulse brings about a set of transition
probabilities to many possible future states (right).

In designing our model, we place high priority on speed and fitting accuracy.
One of the beginning assumptions is that in every possible device state, the
current can be represented by a linear mixture of two fixed polynomials in
$U$. These two polynomials, which are each estimated from a fit to measurement
data, can be thought of as limiting cases for the highest possible high
resistance state, $I_{\text{HHRS}}(U)$, and lowest possible low resistance
state, $I_{\text{LLRS}}(U)$. The device current in all possible resistance
states is then given by

$I(r,U)=rI_{\text{HHRS}}(U)+(1-r)I_{\text{LLRS}}(U),$ (1)

conveniently reducing the description of the conduction in the material to a
single state variable $0<r<1$. This set of functions are able to be
efficiently evaluated by Horner’s algorithm and serve as a close enough
approximation to the true non-linear conduction behavior for our purposes.

In ReRAM, the overall resistance state as well as the transition behavior is
affected by a vast number of different possible configurations of ionic
defects in the material, giving rise to the observed stochastic behavior and
history dependence (Fig. 2). Rather than attempting to describe the ionic
transport physically, we turn instead to measurement data to directly provide
the necessary statistical information. A discrete multivariate stochastic
process based on a Structural Vector Autoregression (SVAR) model is fit to the
data and used to generate latent variables that guide the state evolution of
simulated memory cells. As a cell is exposed to voltage signals, new terms of
the SVAR model are realized by a sum of easily computable linear
transformations of past states and pseudorandom vectors.

As an overview, the experimental and simulation approach that will be
elaborated in this section can be shortly summarized as follows:

  1. 1.

A fabricated ReRAM cell is experimentally driven through a large number of
resistance cycles by applying a continuous periodic voltage signal while
measuring the resulting current.

  2. 2.

A time series of feature vectors, $\bm{x}_{n}$, composed of resistance values
and switching threshold voltages, is extracted from each of the measured
cycles.

  3. 3.

A discrete stochastic process, $\bm{x}^{*}_{n}$, is constructed to enable
generation of simulated feature vectors that reproduce the measured
distributions as well as the long range correlation structure of $\bm{x}_{n}$.

  4. 4.

An array of simulated cells are instantiated according to independent
realizations of $\bm{x}^{*}_{n}$ to represent cycle-to-cycle variations,
together with a random scaling vector $\bm{s}_{m}$ to represent device-to-
device variations.

  5. 5.

Two programming methods are exposed for each cell; one to apply voltages and
another to make realistic current readouts. Applied voltages above the
generated thresholds alter the device state, following an empirical structure
which encodes the resistance transition behavior and allows access to a range
of resistance states. Each voltage driven resistance cycle triggers the
generation of new stochastic terms from $\bm{x}^{*}_{n}$, which govern the
progression to future states.

###  4.1 Data collection

For the purposes of stochastic modelling, electrical measurement data is
needed that capture relevant information about the internal state of a memory
cell and its variation cycle-to-cycle (CtC) and device-to-device (DtD).
However, ReRAM measurements performed at operational speed typically make
exclusive use of rectangular voltage pulse sequences, which yield very little
useful state information. On the other hand, measurements applying
continuously swept voltage signals while sampling the resulting current are
more suitable because much more information is collected each cycle, such as
switching threshold voltages, current-voltage nonlinearity, resistance states,
and transition behavior.

Conventionally, measurements employing voltage sweeps are carried out using
the source measure units (SMUs) of commercial semiconductor parameter
analyzers (SPAs). However, SMUs make heavy use of averaging to measure noisy
signals at high resolution and thus sample too slowly to collect cycling data
in a meaningful quantity. Furthermore, because two-terminal switching devices
are prone to electrical instability and runaway transitions, voltage sweeping
measurements usually require integrated current limiting transistors to avoid
destruction or rapid degradation of the cell. This presents a significant
fabrication overhead and limits the materials available for study. In light of
these challenges, the input data for the present stochastic model was acquired
using a custom measurement technique, introduced in detail in a recent
publication [30]. The setup uses an external current-limiting amplifier
circuit to allow for collection of sweeping measurements at over six orders of
magnitude higher speeds than SMUs, while also eliminating the cumbersome
requirement of on-chip current limiting.

The ReRAM cell used for measurement of cycling statistics (Fig. 3) was
fabricated using a combination of atomic layer deposition (ALD) and physical
vapor deposition (PVD) [31]. The material stack from bottom to top was TiN /
ALD HfO2 (10 nm) / PVD Ti (20 nm) / PVD TiN (100 nm). The device was
electrically isolated with contact pads leading directly to the top and bottom
device electrodes, using no access transistor or added series resistance.
Using a fixed 100 $\upmu$A current limit in the SET polarity, the pristine
cell was electroformed by application of 100 $\upmu$s duration triangular
pulses with incrementally increasing amplitude until a current jump was
recorded near 3 V. For all subsequent cycling, a 1.5 V amplitude 10 KHz
triangular waveform was applied. The cell was first exercised for $2.4\times
10^{6}$ cycles before $10^{6}$ additional cycles were collected for analysis.
Current ($I$) and voltage ($U$) waveforms were simultaneously recorded with
8-bit resolution and with a sample rate of 1,042 samples per cycle.

Figure 3: Scanning electron micrographs of the ReRAM cell design used for
electrical measurement. (A) shows a cross-section of the cell, and (B) shows a
zoom-in of the resistive memory between metalization layers M4 and M5.

The measured current array was smoothed with a moving average filter to
improve the quality of the raw data before further analysis. An adaptive
rectangular window size was used in order to preserve current steps in the
signal, with the maximum window size of 25 samples gradually reducing to a
minimum of 3 samples at the pre-detected locations of SET transitions of each
cycle. After smoothing, the contiguous $I$ and $U$ waveforms were split into
indexible cycles at most positive value of the periodic applied voltage (see
Fig. 4).

Figure 4:  Measured time dependence of $I$ and $U$ waveforms resulting from
the ReRAM cycling experiment. The waveforms are divided into $10^{6}$ indexed
cycles, the first three of which are shown. From this dataset, the periodic
temporal sequence of the states and events of each cycle (HRSn, SETn, LRSn,
RESETn) is extracted and subject to statistical modelling.

Each cycle exhibits the following temporal sequence of states and events: a
high resistance state (HRS), a transition (SET) out of the HRS into the
following low resistance state (LRS), and finally another transition (RESET)
into the next HRS. Current vs. voltage ($I,U$) plots for a subset of the
collected cycles are shown in Fig. 5, which highlights the significant
stochastic CtC variations. The observed characteristics are typical for ReRAM
subjected to voltage controlled sweeps — on average, there is relatively
higher voltage non-linearity in the HRS than in the LRS, and the SET
transitions are abrupt with respect to the applied voltage, while the RESET
transitions proceed relatively gradually over a voltage range of approximately
700 mV.

Figure 5: A subset of the $10^{6}$ measured ($I,U$) cycles used as input to
the stochastic model. The black arrowed path shows the average ($I,U$) curve
and its temporal direction. Different cycle indices are represented by colored
paths, which show significant statistical variation.

###  4.2 Feature extraction

The full $I,U$ cycling measurement just described consists of over 16 GB of
numerical data and would not be practical to model on a point-by-point basis.
Therefore, we aim to compress the dataset while retaining enough information
such that the full ($I,U$) characteristics can be approximately reconstructed
from the compressed representation. Accordingly, the full dataset is reduced
to a vector time series of distinguishing features of each cycle. Four scalar
features were chosen for extraction: the value of the HRS, $R_{H}[\Omega]$,
the SET threshold voltage, $U_{S}[V]$, the value of the LRS, $R_{L}[\Omega]$,
and the RESET voltage, $U_{R}[V]$. We denote the series as

$\bm{x}_{n}=\begin{bmatrix}R_{H,n}\\\ U_{S,n}\\\ R_{L,n}\\\
U_{R,n}\end{bmatrix}=\begin{bmatrix}R_{H}\\\ U_{S}\\\ R_{L}\\\
U_{R}\end{bmatrix}_{n},$ (2)

where $n=\\{1,2,\ldots,10^{6}\\}$ is the set of cycle indices. The feature
vector elements, whose precise definition follows, are chronologically ordered
from top to bottom as they occur in the measurement dataset.

The SET voltage $U_{S}$, or the voltage where the cell resistance abruptly
decreases, is extracted from each cycle as the absolute value of the linearly
interpolated $U$ corresponding to the first level crossing of
$I=-50~{}\upmu\text{A}$. The RESET voltage $U_{R}$, defined as the voltage
where the reset process begins, is determined from the $I$ datapoints by peak
detection using simple comparison of neighboring samples. Here, only the
increasing section of the voltage sweep with $U>0$ is considered. The voltage
corresponding to the first encountered peak with prominence $\geq 5~{}\upmu$A
is taken as the RESET voltage. If no peak satisfies this criterion, the peak
with maximum prominence is taken instead.

The device current for any static state is approximated in our model as a
polynomial function of the applied voltage. The values of $R_{H}$ and $R_{L}$
are likewise extracted from least squares polynomial fits to appropriate
subsets of the measured $(I,U)$ data of each cycle. The HRS is fit with a 5th
degree polynomial on the decreasing $U$ sweep in the variable range
$U_{S}+0.1~{}\text{V}~{}\leq~{}U~{}\leq~{}1.5~{}\text{V}$ and
$-25~{}\upmu\text{A}~{}\leq~{}I~{}\leq~{}80~{}\upmu\text{A}$, and the LRS is
fit with a 3rd degree polynomial on the increasing part of the $V$ sweep in
the range $-0.7~{}\text{V}~{}\leq~{}U~{}\leq~{}U_{R}-0.05~{}\text{V}$ and
$-80~{}\upmu\text{A}~{}\leq~{}I~{}\leq~{}120~{}\upmu\text{A}$. The fits are
constrained such that the 0th order coefficient equals 0 A, and the 1st order
coefficient is $\geq$1 nA/V. The values of $R_{H}$ and $R_{L}$ are then
defined as the static resistance of the respective polynomials at a fixed
voltage $U_{0}=200$ mV.

An overview of the result of this feature extraction is given in Fig. 6. The
$10^{6}$ cycles proceeded without significant long-term drift from the overall
mean value,

$\bar{\bm{x}}_{n}=\begin{bmatrix}166.5~{}\text{k}\Omega\\\ 0.85~{}\text{V}\\\
8.2~{}\text{k}\Omega\\\ 0.72~{}\text{V}\end{bmatrix},$ (3)

but with significant variations in each feature between cycles. A prominent
characteristic of this data is that it is strongly correlated over long cycle
ranges, as quantified in Fig. 14. The asymmetric marginal distributions for
each of the features were very well resolved due to the large number of
samples, and they did not accurately converge to any analytical probability
density function (PDF) in common use, including the normal and log-normal.

Figure 6: A view of the feature vector time series extracted from each of
$10^{6}$ measured ($I,U$) cycles. Each feature, which represents either a
resistance state or a switching voltage, has its marginal histogram shown on
the right.

###  4.3 Stochastic Modelling

This section will introduce the statistical methods used to model the internal
states of an array of synaptic ReRAM devices, including CtC and DtD
variability effects. The handling of voltages applied to the cells as well as
the simulation of realistic readouts of the resistance states will also be
established. To help orient the reader, the overall structure of the
generative model that will be described is provided in advance in Fig. 7.

$a$ $\bm{\Sigma}$ $\bm{\gamma}$ $\bm{A}$ $\bm{B}$ $\bm{C}_{i}$ $U_{\text{max}}$ $\bm{p}$ $U_{\text{read}}$ $\Delta f$ $I_{\text{min}}$ $I_{\text{max}}$ $n_{\text{bits}}$ $\widehat{\bm{s}}_{m}$ $\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}$ $\bm{s}_{m}$ $\bm{\epsilon}_{n}$ $\widehat{\bm{x}}^{*}_{n}$ $\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}$ $\bm{x}^{*}_{n}$ $\bm{y}^{*}_{m,n}$ $r$ $U_{a}$ $I_{\text{read}}$ $\sigma_{I}$ $I_{\text{ADC}}$ $\forall m\in\\{1,\ldots,M\\}$$\forall n\in\\{1,\ldots,N\\}$$n\mapsto\\{n+1,\ldots,n+p\\}$ $a$ |  Scale factor for $\bm{\Sigma}$ | $\widehat{\bm{s}}_{m}$ |  $\mathcal{N}(\bm{0},a\Sigma)$ sample
---|---|---|---
$\bm{\Sigma}$ | Device covariance | $\bm{\Gamma}^{\raisebox{0.48221pt}{$\scriptscriptstyle\text{-}1$}}$ | Inverse normalizing map
$\bm{\gamma}$ |  Coefficients for $\bm{\Gamma}^{\raisebox{0.48221pt}{$\scriptscriptstyle\text{-}1$}}$ | $\bm{s}_{m}$ | Device scale vector
$\bm{A}$ | SVAR contemp. parameters | $\bm{\epsilon}_{n}$ | White noise vector
$\bm{B}$ | SVAR noise amplitudes | $\widehat{\bm{x}}^{*}_{n}$ | SVAR process
$\bm{C}_{i}$ | SVAR lag parameters | $\bm{x}^{*}_{n}$ | Median cycling process
$U_{\text{max}}$ | Max. voltage applied | $\bm{y}^{*}_{m,n}$ | Device cycling process
$\bm{p}$ | HHRS, LLRS coefficients | $r$ | Device state variable
$U_{\text{read}}$ | Readout voltage | $U_{a}$ | Applied voltage pulse(s)
$\Delta f$ | Readout bandwidth | $I_{\text{read}}$ | Readout current (noiseless)
$I_{\text{min}}$ | Min. current of ADC | $\sigma_{I}$ | Readout noise amplitude
$I_{\text{max}}$ | Max. current of ADC | $I_{\text{ADC}}$ | ADC readout
$n_{\text{bits}}$ | Number of ADC bits |  |
Constant node Stochastic node Deterministic node Observed node Function Plate
0D 1D 2D 3D  Figure 7:  Graphical model depicting the relationships between
all parameters and latent variables involved in the stochastic synapse model.
Plate notation is used to represent $N$ switching cycles of $M$ devices, each
yielding an observed readout current. The dotted recurrent arrow denotes a
connection to each of the $p$ following frames, as needed by the history
dependent stochastic process.

####  4.3.1 Cycle-to-cycle (CtC) variations

In seeking to represent the input time series $\bm{x}_{n}$ with a stochastic
process, the main goals are to recreate the marginal distributions as well as
the correlation structure of its vector components. To achieve the first goal
with high generality, we use an approach based on transformation of the
measured densities to and from the standard normal distribution
$\mathcal{N}(0,1)$. This way, a single process can be used to achieve any set
of marginals presented by the input data, with the relatively unrestrictive
requirement that this base process generates normal marginals. Notationally,
we define and apply an invertible, smooth mapping
$\bm{\Gamma}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ that normalizes the
marginal distributions of the vector components,

$\bm{x}_{n}=\begin{bmatrix}R_{H}\\\ U_{S}\\\ R_{L}\\\
U_{R}\end{bmatrix}_{n}\xrightarrow{\bm{\Gamma}}\begin{bmatrix}\widehat{R}_{H}\\\
\widehat{U}_{S}\\\ \widehat{R}_{L}\\\
\widehat{U}_{R}\end{bmatrix}_{n}=\widehat{\bm{x}}_{n},$ (4)

where a hatted variable signifies that it is distributed as
$\mathcal{N}(0,1)$. We then construct a base process
$\widehat{\bm{x}}^{*}_{n}$ whose marginals are normal, and finally transform
its output back to the original data distributions via the inverse map
$\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}$. The
overall process $\bm{x}^{*}_{n}$ is thus defined,

$\widehat{\bm{x}}^{*}_{n}=\begin{bmatrix}\widehat{R}_{H}^{*}\\\
\widehat{U}_{S}^{*}\\\ \widehat{R}_{L}^{*}\\\
\widehat{U}_{R}^{*}\end{bmatrix}_{n}\xrightarrow{\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}}\begin{bmatrix}R_{H}^{*}\\\
U_{S}^{*}\\\ R_{L}^{*}\\\ U_{R}^{*}\end{bmatrix}_{n}=\bm{x}^{*}_{n},$ (5)

where a star indicates a generated random variable to distinguish from
variables originating from measurement data.

This type of density transformation procedure is a widely used technique for
working with arbitrary distributions, which finds application in a variety of
fields and can be constructed in many different ways [32, 33]. While the
transformation is trivially constructed in the case where the target quantile
function and its inverse are each analytically defined, we do not make this
assumption in the present scenario. A simple numerical method in this case is
a so-called quantile transform, where the input and output quantile functions
are each discretely sampled and the transformation is defined through a direct
map between bins or through interpolation. The main requirement for
$\bm{\Gamma}$ in our model, however, is that its inverse (Eq. 5) is easy to
evaluate without causing cache misses due to memory access, thus it is
preferable to avoid referencing and interpolation of large look-up tables. The
forward transformation (Eq. 4), on the other hand, only needs to be computed
once for model fitting and is not used for the generating process. We
therefore define
$\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}$ as
essentially a quantile transform, operating on each feature independently,
that is evaluated from a fit of the quantiles to a specific analytic function.
Namely,

$\bm{\bm{\Gamma}}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\widehat{\bm{x}}_{n})=\exp\begin{bmatrix}[l]\gamma_{1}(\widehat{R}_{H,n})\\\
\gamma_{2}(\widehat{U}_{S,n})\\\ \gamma_{3}(\widehat{R}_{L,n})\\\
\gamma_{4}(\widehat{U}_{R,n})\end{bmatrix}=\bm{x}_{n},$ (6)

where $\gamma_{1}$-$\gamma_{4}$ are each 5th degree polynomials, and the
exponential function is applied element-wise. The coefficients of the
polynomials are fit to standard normal quantiles vs. those of the respective
(log) features, sampled at 500 equally spaced values between 0.01 and 0.99.
The fitted polynomials are checked for monotonicity within four standard
deviations above and below zero, and the forward transformation,

$\bm{\bm{\Gamma}}(\bm{x}_{n})=\begin{bmatrix}[l]\gamma_{1}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\log{R_{H,n}})\\\
\gamma_{2}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\log{U_{S,n}})\\\
\gamma_{3}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\log{R_{L,n}})\\\
\gamma_{4}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\log{U_{R,n}})\end{bmatrix}=\widehat{\bm{x}}_{n},$
(7)

is computed using numerical inverse of the $\gamma$ polynomials. A
visualization of the function $\bm{\Gamma}$ as well as the marginal histograms
corresponding to input series $\bm{x}_{n}$ and output series
$\widehat{\bm{x}}_{n}$, are shown in Fig. 8.

Figure 8: Visualization of the invertible normalizing transformation
$\bm{\Gamma}$ that is applied to the measured feature vectors before fitting
with a base stochastic process. The left column shows the marginal PDFs of the
vector time series $\bm{x}_{n}$ extracted from measurement. The center column
shows the input and output quantile-quantile plots with the fitted log-
polynomial function used to transform the distributions (here, $Q$ denotes the
quantile function of its argument). The right column is the result of applying
$\bm{\Gamma}$ to the input data, producing $\widehat{\bm{x}}_{n}$ whose
elements are normally distributed.

Now that we have transformed the input measurement data into a normalized
vector time series $\widehat{\bm{x}}_{n}$, a suitable stochastic process will
be chosen for fitting. This process should serve as a useful approximation to
the true physical mechanisms that generated the data, capturing the long-range
correlation structure of the observed features. Time series analysis is
broadly used across scientific and engineering domains, but despite its
applicability to the rich statistical behavior displayed by resistive
switching devices, device models have not yet widely employed dependent
stochastic processes. Many models and analyses assume for convenience that
features are independently and identically distributed according to a normal
or lognormal PDF [20, 34]. However, there is not a strong theoretical basis
for this assumption in a highly nonlinear and path-dependent system based on
continuous evolution of conducting filaments. Dependent stochastic processes,
on the other hand, more appropriately allow for a description of the
dependence of future states on past states.

Simple models in the category of Markov chains have been considered as
generating processes for memory cells. A rudimentary example is a
1-dimensional random walk process, where each future state is computed as a
random additive perturbation on the previous state [22]. While random walk
represents a reasonable short range approximation, it has the well known
property that the expected absolute distance between the initial value and the
$N$th value is proportional to $\sqrt{N}$ for large $N$, causing the process
to eventually drift to unphysical values without the use of artificial
constraints.

Autoregressive (AR) models are simple univariate processes sharing some
characteristics of random walk, but based additionally on a deterministic
linear dependence on past observations. Each new term of an AR($p$) (AR of
order $p$) model is computed by linear combinations of $p$ previous (lagged)
values together with a noise term, producing processes that are wide-sense
stationary and mean-reverting within suitable parameter ranges [35, 36]. The
few times they have appeared in the literature, low order models like AR(1)
and AR(2) were used to describe state variables independently (e.g. a sequence
of high and/or low resistance states) [37, 38]. Here we pursue a more
comprehensive statistical description of the interrelations between the
different variables contained in the vectors $\widehat{\bm{x}}_{n}$ which
takes into account long range correlations $p\gg 1$. This is enabled by the
use of a VAR($p$) model (vector AR of order $p$), which is the multivariate
counterpart of the AR model applicable to discrete vector time series [35,
36].

We adopt in particular a Structural VAR (SVAR) formulation of the model, which
is a factorization that makes the relationships between the contemporaneous
(same index) variables explicit. The model has the form

$\bm{A}\widehat{\bm{x}}^{*}_{n}=\sum_{i=1}^{p}\bm{C}_{i}\widehat{\bm{x}}^{*}_{n-i}+\bm{B}\bm{\epsilon}_{n},$
(8)

where $\bm{A}$, $\bm{B}$, and $\bm{C}_{i}$ are $4\times 4$ matrices of model
parameters, and $\bm{\epsilon}_{n}$ is a 4-dimensional standard white noise
process. With this formulation we impose a general structure of causal
ordering for the generated random variables consistent with the chronological
chain of measurement events. Within this structure, each variable may have a
causal and deterministic effect on all future variables within range $p$, as
visualized by the graph of Fig. 9. The size of these effects are all subject
to fitting via the coefficients of the model. Constraints on the structural
parameters,

$\bm{A}=\begin{bmatrix}1&0&0&0\\\ A_{21}&1&0&0\\\ A_{31}&A_{32}&1&0\\\
A_{41}&A_{42}&A_{43}&1\end{bmatrix},\bm{B}=\begin{bmatrix}B_{11}&0&0&0\\\
0&B_{22}&0&0\\\ 0&0&B_{33}&0\\\ 0&0&0&B_{44}\end{bmatrix}$ (9)

enforce the desired causal structure while assuming an uncorrelated noise
driving process. Model fitting was performed using the Python statsmodels
package [39], wherein a VAR($p$) model is first fit by ordinary least squares
regression, and a maximum likelihood estimate is then used to determine the
structural decomposition.

$\widehat{R}^{*}_{H,n-1}$ $\widehat{U}^{*}_{S,n-1}$ $\widehat{R}^{*}_{L,n-1}$
$\widehat{U}^{*}_{R,n-1}$ $\widehat{R}^{*}_{H,n}$ $\widehat{U}^{*}_{S,n}$
$\widehat{R}^{*}_{L,n}$ $\widehat{U}^{*}_{R,n}$ $\epsilon_{1,n}$
$\epsilon_{2,n}$ $\epsilon_{3,n}$ $\epsilon_{4,n}$
-0.0020.0370.0210.0430.9840.111-0.0110.0280.0570.945-0.1390.0230.0100.1530.9080.1800.070-0.0080.0850.921$\widehat{\bm{x}}^{*}_{n-1}$$\widehat{\bm{x}}^{*}_{n}$
Figure 9: A weighted graph displaying the causal structure of the utilized
SVAR($p$) process, showing the nearest temporal contributions to realizations
of the random vector $\widehat{\bm{x}}^{*}_{n}$. Arrow weights show the model
parameters contained in $\bm{A}$, $\bm{B}$ and the upper triangular part of
$\bm{C}_{1}$ when fit with $p=100$. The actual SVAR($p$) model uses many more
connections than shown ($16p+10$), so that each variable is impacted by all
past values of all other variables within cycle range $p$.

####  4.3.2 Device-to-device (DtD) variations

So far, we have only considered the statistical modelling of the cycling
process of a single memory cell. However, the purpose of the presented model
is to simultaneously simulate a large number of cells in a network. Individual
memory devices on a wafer generally show statistical variations, mainly
arising due to defects and non-uniformities in fabrication [40, 41]. These DtD
variations depend strongly on the particular lithography processes and
materials used. They can also originate from intrinsic factors and are
influenced by conditions during the electroforming of each cell [42, 43].
Because of the potential positive or negative impact on network performance,
it is important for the model to account for the DtD variability [41, 44].

The electrical effect of device variability is modelled with each cell using a
modification of the same underlying SVAR cycling process. Device-specific
processes are defined as members of a parametric family of processes, all
based on element-wise scaling of $\bm{x}^{*}_{n}$, where the scaling factors
are themselves random vectors. The specific process is denoted

$\bm{y}_{m,n}^{*}=\bm{s}_{m}\odot\bm{x}_{n}^{*},$ (10)

where $m=\\{1,2,\ldots,M\\}$ is the device index, $\odot$ is the Hadamard
(element-wise) product, and $\bm{s}_{m}$ are $4\times 1$ random vectors drawn
from a fixed distribution at cell initialization.

The distribution of $\bm{s}_{m}$ is chosen so that the features of the median
cycles of different devices are distributed and correlated in the same way as
the measured cycling data $\bm{x}_{n}$. This choice reflects that the
covariations of switching features DtD arise in the same physical system with
causes and effects that are comparable to those of the CtC variations. To this
end, random vectors $\widehat{\bm{s}}_{m}$ are drawn from a multivariate
normal (MVN) distribution and
$\Gamma^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}$ is then reused
to map them to the measured CtC distribution,

$\bm{s}_{m}=\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\widehat{\bm{s}}_{m})\oslash\bm{\Gamma}^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}(\bm{0}),\text{
where }\widehat{\bm{s}}_{m}\sim\mathcal{N}(\bm{0},a\bm{\Sigma}).$ (11)

Here, the denominator of the Hadamard division ($\oslash$) sets the median
scale vector to the identity, ${\bm{\Sigma}=\text{cov}(\widehat{\bm{x}}_{n})}$
is the sample covariance of the normalized measurement data, and $a$ is a free
scalar parameter providing adaptability to different DtD covariance levels. A
robust determination of $a$ requires measurement of many switching cycles
across a large number of devices of interest. Values in the range
$a\in[1,1.5]$ approximately correspond to published DtD measurement samples
[40, 41], but improved processing and electroforming procedures may justify
the use of $a<1$.

####  4.3.3 Control logic

As components of a network, each simulated cell possesses a resistance state
that encodes the weight of a connection. Voltage pulses directly applied to
the cells are used to produce resistance state transitions to update the
weights. In this model, applied voltage pulses are distinguished only by a
scalar amplitude $U_{a}$, whether they are in fact square waveforms or they
have a more complex shape of an action potential. Although ReRAMs are known to
be highly time-dependent devices [45], we assume here that the duration of the
pulses is appropriately matched to the experimental timescale, such that a
simulated voltage pulse of a given amplitude produces an effect comparable to
the experimental voltage sweep at the instant it reaches that same amplitude.
Possible state modifications in response to an input pulse is computed with
respect to $I,U$ sweeps that are reconstructed from each stochastic feature
vector generated for each cycle as illustrated in Fig 10.

Figure 10:  Conduction polynomials and threshold voltages allow reconstruction
of ($I,U$) cycles from generated feature vectors. Simulated resistance
switching is such that the conduction state $I(r,U)$ induced by an applied
voltage $U_{a}$ intersects the reconstructed cycle at $U=U_{a}$. For visual
simplicity, the cycle shown begins and ends in the same HRS
($R_{H,n}=R_{H,n+1}$).

As previously specified in eq. 1, every possible electrical state of a device
is assumed to correspond to a polynomial $I(U)$ dependence parameterized by a
state variable $r$. It is straightforward to calculate that the state variable
for a curve passing through an arbitrary ($I,U$) point is uniquely given by
the function

$r(I,U)=\frac{I_{\text{LLRS}}(U)-I}{I_{\text{LLRS}}(U)-I_{\text{HHRS}}(U)}.$
(12)

Therefore the state variable corresponding to any static resistance level $R$
(evaluated at $U_{0}$) can be calculated using

$r(R)=\frac{I_{\text{LLRS}}(U_{0})-U_{0}R^{\raisebox{0.60275pt}{$\scriptscriptstyle\text{-}1$}}}{I_{\text{LLRS}}(U_{0})-I_{\text{HHRS}}(U_{0})}.$
(13)

The $I(U)$ curves for the electrical states corresponding to each cycle’s HRS
and LRS, hereafter called $I_{\text{HRS,n}}(U)$ and $I_{\text{LRS,n}}(U)$, are
defined according to equations (1) and (13) such that their static resistance
equals the respective value of $R^{*}_{H,n}$ and $R^{*}_{L,n}$.

Transitions between HRS and LRS states in response to an applied pulse
amplitude $U_{a}$ follow an empirically motivated structure, represented by
the flow chart of Fig. 11. The SET transition for the $n$th cycle HRSn
$\rightarrow$ LRSn may occur for negative voltage polarities and follows a
simple threshold behavior, fully and instantaneously transitioning the first
time a voltage pulse with amplitude $U_{a}\leq U^{*}_{S,n}$ is applied. In
contrast, the RESET transition LRSn $\rightarrow$ HRSn+1 occurs gradually in
the positive polarity with increasing $U_{a}$ in the range
$U^{*}_{R,n}<U_{a}\leq U_{\text{max}}$, where $U_{\text{max}}=1.5$ V is the
maximum voltage applied in the voltage sweeping measurement. A transition
curve $I_{\text{RESET,n}}(U)$ is defined to connect the ($I,U$) points of the
two limiting states where the RESET transition begins and ends. The functional
form of the transition curve is chosen to be the parabola with boundary
conditions

$\displaystyle I_{\text{RESET,n}}(U^{*}_{R,n})$ $\displaystyle=$
$\displaystyle I_{\text{LRS,n}}(U^{*}_{R,n})$ (14) $\displaystyle
I_{\text{RESET,n}}(U_{\text{max}})$ $\displaystyle=$ $\displaystyle
I_{\text{HRS,n+1}}(U_{\text{max}})$ (15)
$\displaystyle\left.\frac{dI_{\text{RESET,n}}}{dU}\right|_{U=U_{\text{max}}}$
$\displaystyle=$ $\displaystyle 0.$ (16)

When a voltage pulse in the RESET range is applied, an intermediate resistance
state (IRS) results which is calculated with reference to the transition curve
such that $I(r,U_{a})=I_{\text{RESET,n}}(U_{a})$. Additional RESET pulses with
larger amplitudes may be applied to incrementally increase the cell
resistance, with HRSn+1 being reached only if $U_{a}\geq U_{\text{max}}$,
after which no further RESET switching is possible for the $n$th cycle. After
either partial or full RESET, the resistance may only decrease again by
entering the following LRSn+1 with a voltage pulse meeting the SET criterion
$U_{a}\leq U^{*}_{S,n+1}$.

$n$ | $\leftarrow$ | $1$
---|---|---
$r$ | $\leftarrow$ | HRS1
$U_{R}$ | $\leftarrow$ | $U^{*}_{R,1}$
Cell Initialization$U_{a}\leftarrow$ Pulse Input$U_{a}>U_{R}$$U_{a}\leq U^{*}_{S,n}$in HRSin LRS$U\geq U_{\text{max}}$in LRS$r\leftarrow$ HRSnFull RESET $r$ | $\leftarrow$ | IRSn($U_{a}$)
---|---|---
$U_{R}$ | $\leftarrow$ | $U_{a}$
Partial RESET$r\leftarrow$ LRSn SET$n\leftarrow n+1$FTTFFTTFFFTT Figure 11:
Logical flow chart showing how applied voltage pulses affect the state of each
cell during simulation. Following the experimental observations, SET processes
always occur abuptly below a threshold voltage, while partial switching is
induced for a range of RESET voltages, with intermediate states bounded for
cycle $n$ by resistance values between $R_{L,n}$ and $R_{H,n+1}$. As
resistance cycling progresses, later terms of the stochastic driving process
are used for limiting resistance states and threshold voltages. Pulse
amplitudes not producing a state change are efficiently disregarded.

####  4.3.4 Readout

Simulated current measurements (readouts) for each individual cell can be
generated given an arbitrary readout voltage input $U_{\text{read}}$. The
noise-free current level simply corresponds to evaluation of
$I(r,U_{\text{read}})$ for each cell. In any real system, however, current
readouts are accompanied by measurement noise, which may impact system
performance and even present a fundamental bottleneck. Furthermore, in digital
systems current readouts are converted to finite resolution by analog to
digital converters (ADCs). Due to constraints of power consumption and chip
area, ADC resolution is often limited such that digitization is the dominant
contributor to the total noise [46]. Many additional noise sources can be
considered, such as $1/f$ noise [47], but at minimum the Johnson-Nyquist noise
and the shot noise should be included because they represent a lower bound of
noise amplitude impacting all systems.

To account for measurement noise, each individual current readout includes an
additive noise contribution drawn from a normal distribution. The noise
amplitude is approximated from the Nyquist and Schottky formulas,

$\sigma_{I}=\sqrt{\frac{4k_{B}TI_{\text{read}}\Delta
f}{U_{\text{read}}}+2qI_{\text{read}}\Delta f},$ (17)

where $\Delta f$ is the noise equivalent bandwidth, $k_{B}$ is the Boltzmann
constant, $T=300$ K is the temperature, $q$ is the electron charge,
$I_{\text{read}}$ is the noiseless current readout, and $U_{\text{read}}$ is
the voltage used for readout. The total current is then ideally digitized with
an adjustable resolution $n_{\text{bits}}$ between adjustable minimum
$I_{\text{min}}$ and maximum $I_{\text{max}}$ current levels.

###  4.4 Program implementation

To facilitate investigations of neuromorphic systems, model implementations
designed to simulate arrays of devices were developed in the Julia programming
language. Julia is a modern high-level language that is focused on performance
and that provides an advanced ML and scientific computing ecosystem. Julia
programs compile to efficient native code for many platforms via the LLVM
compiler infrastructure, and a cursory analysis indicated that single threaded
CPU performance of a Julia implementation is up to 5,000 times faster than a
Python implementation. Furthermore, as modern computational resources are
highly parallel, Julia’s support for CPU multi-threading and GPU programming
through CUDA.jl [48] is an important advantage.

All model parameters corresponding to the device characterized in this
article, including different possible SVAR model orders, $p\in[1,200]$, are
stored in a binary file which is read in by the program at startup. Each
instantiated cell stores state information and $p$ cycles of history using
primarily 32-bit floating point numbers. The total memory footprint grows
linearly with the chosen model order and is approximately $16p+56$ bytes per
cell. A reduced form VAR process is used to compute realizations of
$\bm{x}^{*}_{n}$, which are lazily evaluated along with the parabolic
transition polynomials if and when they are needed. The majority of the
necessary runtime computations are formulated as matrix multiplications, which
are heavily optimized operations across many different contexts.

The present release contains two model implementations in order to suit a wide
variety of computing platforms and use cases [49]. The first is a CPU
optimized version wherein the cells of an array are individually addressable
for read/write operations. These operations are naturally parallelized for
multi-core processors by partitioning the cells and assigning each partition
to independent threads of execution. The second implementation is a GPU
accelerated version compatible with CUDA capable GPUs. This version uses a
vectorized data structure and parallel array abstractions to take advantage of
the implicit parallelism programming model of CUDA.jl. Here, all defined cells
are always accessed simultaneously, with each read/write operation employing
optimized linear algebra GPU kernels. While the GPU implementation integrates
well with other ML components risiding in GPU shared memory and achieves
higher throughput per cell for large parallel operations, the CPU
implementation obtains higher update rates for sparse operations commonly
encountered in large-scale models [50, 51].

##  5 Results

As shown visually in the scatterplot of Fig. 12, the stochastic process
$\bm{x}^{*}_{n}$ generates data that closely resemble the measurement data
$\bm{x}_{n}$. To quantify the difference between the generated distributions
and the empirical distributions, the first Wasserstein distance [52] was
calculated element-wise and averaged across 100 realizations of
$\bm{x}_{n}^{*}$ with length $10^{6}$. The result,

$\mkern 1.5mu\overline{\mkern-1.5mu\bm{W}\mkern-1.5mu}\mkern
1.5mu_{1}(\bm{x}_{n},\bm{x}_{n}^{*})=\begin{bmatrix}5146~{}\Omega\\\
937~{}\upmu\text{V}\\\ 20~{}\Omega\\\ 356~{}\upmu\text{V}\end{bmatrix},$ (18)

is much smaller than the mean feature vector, $\bar{\bm{x}}_{n}$ (Eq. 3), and
independent of the chosen model order. This shows that the goal of reproducing
the measurement distributions is well achieved for the input dataset by using
the described method of probability density transformation.

Figure 12: Comparison of feature time series extracted from measurement data
and those generated by the SVAR-based model. The compared features converge to
effectively equivalent distributions and the short-range behavior is
qualitatively similar across thousands of cycles.

Simulations of full ($I,U$) cycling measurements [Fig. 13(A)] show close
similarity with the measurement data of Fig. 5. Multi-resistance-level
capability is also demonstrated by a similar simulation involving partial
RESET operations by changing the maximum voltage applied [Fig. 13(B)]. The
dependence of the resulting HRS value on the applied voltage reproduces a non-
linear characteristic comparable to experimental findings [53, 54].

Figure 13:  Two example simulations involving repeated cycling of a single
device. Voltage pulse sequences were applied with varying amplitude following
a triangular envelope, and the ($I,U$) characteristic of each cycle is plotted
in a different color. Subplot (A) shows 300 consecutive cycles between the
full voltage range $\pm$1.5 V, with a readout performed after every pulse
(inset). Subplot (B) demonstrates multilevel capability with 300 cycles
between -1.5 V and maximum voltage that increases each cycle, from 0.7 V to
1.5 V. Readouts following each cycle are shown in the inset. In each case,
readouts were simulated using a fixed $U_{\text{read}}=$ 200 mV, including
noise and 4-bit quantization between $I_{\text{min}}=0~{}\upmu$A and
$I_{\text{max}}=40~{}\upmu$A.

While a full structural analysis of the fitted SVAR($p$) model parameters
($\bm{A},\bm{B},\bm{C}_{i}$) will not be presented here, a few aspects are
worthy of note. For the fit corresponding to the particular device and
measurement described here, the white noise terms are by far the dominant
contributors to all four modelled features. The contemporaneous terms
($\bm{A}$) and first order ($\bm{C}_{1}$) terms are the next most significant,
which indicates that the most recent cell history is most relevant for
generating the proceeding states. Nevertheless, input data correlations
persist for many cycles, and the generating process $\bm{x}^{*}_{n}$
successfully reproduces the overall correlation structure of the data up to at
least $p$ cycle lags, as shown in detail in Fig. 14.

Although no physical effects were explicitly put into the model definition, it
is important to recognize that they are quantitatively captured and put into a
useful statistical context by the SVAR model fitting procedure. As seen in the
graph of Fig. 9, the strongest deterministic coefficients in the fitted model
correspond to the relationships

$\displaystyle\widehat{R}^{*}_{H,i}$
$\displaystyle\xrightarrow{\phantom{-}0.111}\widehat{U}^{*}_{S,i},$ (19)
$\displaystyle\widehat{U}^{*}_{S,i}$
$\displaystyle\xrightarrow{-0.139}\widehat{R}^{*}_{L,i},$ (20)
$\displaystyle\widehat{R}^{*}_{L,i-1}$
$\displaystyle\xrightarrow{\phantom{-}0.153}\widehat{R}^{*}_{L,i},$ (21)
$\displaystyle\widehat{R}^{*}_{L,i}$
$\displaystyle\xrightarrow{\phantom{-}0.180}\widehat{U}^{*}_{R,i}.$ (22)

Comparable relationships between switching variables have been identified and
discussed in physics-based models and simulations as well as in experimental
studies involving various materials [55, 56, 57, 58, 59, 60, 61]. According to
relation 19, larger starting HRS values tend to contribute to a higher SET
voltage, which is a well known effect due to a reduced driving force for ionic
motion at a given applied voltage, as a larger HRS gives both reduced power
dissipation as well as a reduced electric field in a thicker insulating gap.
The subsequent LRS is strongly affected by the SET voltage (relation 20). This
can be attributed to the runaway nature of the SET transition and a higher
voltage initial condition, and is also connected with the dynamics of the
current limiting circuitry [30]. The LRS value is also strongly correlated
with the value of the previous LRS (relation 21), because of the influence of
the residual filamentary structure from the previous cycle [62]. Lastly,
relation 22 indicates that higher LRS values tend to have larger reset
voltages, which has to do with a balance of factors influencing filament
dissolution, including temperature and drift. This balance depends on the cell
materials, operating regime, and internal series resistance [63].

Figure 14:  Auto- and cross-correlations of the normalized feature vector
components, showing the pearson coefficients $\rho_{X,Y}$ of the variables
specified in the subplot columns $X$ and rows $Y$ as a function of lag $l$.
Row variables are lagged with respect to column variables, as denoted by the
lag operators $L_{X}$. A comparison between measurement data and data
generated from SVAR(30) shows extremely close agreement up to cycle range 30.
For lags larger than the chosen model order, some of the correlations of
$\widehat{\bm{x}}^{*}$ decay more quickly than $\widehat{\bm{x}}$.

###  5.1 Benchmarks

As a benchmark of the throughput of write operations, repeated resistance
cycling was induced on arrays of simulated cells under varying conditions. In
each case, voltage pulse sequences to be applied to all defined cells were
generated prior to the benchmarks, consisting of amplitudes $\pm$1.5 V with
alternating polarity. Defined as such, every pulse drives each cell through a
transition into its next HRS or LRS. The read operation was benchmarked
separately under equivalent conditions, reading out the entire array using a
fixed readout voltage of $U_{\text{read}}=0.2~{}V$.

The CPU benchmark was performed using an Intel Xeon Silver 4116 CPU, varying
the cell array size $M$, the order of the VAR process $p$, as well as the
number of threads used to perform the operations in parallel. The resulting
read/write throughputs are summarized in Fig. 15. Write throughputs up to
$2\times 10^{8}$ operations per second (OPS) were obtained, which is
equivalent to 5 ns per individual write operation. Read operations were nearly
an order of magnitude faster than writes, with up to $10^{9}$ OPS or 1 ns per
read operation. Due to the size of necessary matrix multiplications,
increasing the VAR order $p$ incurs a cost of write throughput, with a $p=100$
model running approximately 4$\times$ slower than one with $p=10$. The read
operation, in contrast, shows a negligible dependence on the VAR order.

Figure 15:  Benchmarks of the read/write operation throughput per cell of the
Julia model implementations. In (A) and (C), an array of $2^{20}$
(${\sim}10^{6}$) cells are simulated on the CPU as a function of number of
parallel threads spawned, and the VAR model order $p$. In (B) and (D), the CPU
(32 threads) and GPU implementations are benchmarked versus the cell array
size $M$, with $p=10$.

The GPU accelerated version was benchmarked in an analogous way, using the
same host machine with an NVIDIA TITAN RTX GPU device. The results are shown
in dependence of the cell array size $M$ in Fig. 15(B,D). The GPU
implementation overtakes the CPU above $M=10^{6}$ parallel operations where
the entire array is updated, and achieves 2$\times$ faster updates and
5$\times$ faster readouts for large arrays with $M>10^{7}$. However, CPU
throughput is applicable to subsets of the array, and may retain an advantage
for sparse operations.

##  6 Discussion

In order to assess the potential of emerging synaptic devices, new lightweight
and accurate device models are needed to constitute the millions/billions of
weights used in modern machine learning (ML) models. Candidate memory cells
such as ReRAM are highly non-linear stochastic devices with complex internal
states and history dependence, all of which needs to be explicitly taken into
account. In this article we introduced an efficient generative model for large
synaptic arrays, which closely reproduces the statistical behavior of real
devices.

Taking advantage of a recently developed electrical measurement technique
[30], we systematically fit the model to a dataset that is dense in relevant
information about the device state evolution. Together with this new kind of
measurement, our modelling approach helps complete a neuromorphic design
feedback loop by defining a programmatic connection from the measured behavior
of a fabricated device under the intended operating conditions directly to
fitted model parameters. Probability density transformation of the underlying
SVAR stochastic process gives the model the power to accurately reproduce
nearly arbitrary distribution shapes and covariance structures across the
switching cycles and across the separate devices. These features enable
evaluation of network performance while automatically adapting to a wide
variety of possible future device designs.

We provide parallelized implementations for both CPU and GPU, where up to 15
million cells per GB of available memory can be simulated at once. Benchmarks
show throughputs above three hundred million weight updates per second, which
exceeds the pixel rate of a 30 frames per second video stream at 4K resolution
(3840$\times$2160 pixels). Realistic current readouts including digitization
and noise were also benchmarked, and are approximately an order of magnitude
faster than weight updates. While speeds can be expected to improve with
future optimizations, these benchmarks give a basis for estimating the scope
of applicability of the model to ML tasks.

The implementation and the general concept of this model are naturally
extendable. Although model parameters were adapted here to a specific
HfO2-based ReRAM device, the method is applicable to a variety of other types
of stochastic memory cells such as PCM, MRAM, etc. Four specific switching
features were chosen in this demonstration to reconstruct ($I,U$) cycling
behavior, but additional switching parameters can also be extracted from
measurements and accommodated within this framework. Ideally informed by
statistical measurement data, different functional forms, transition
behaviors, time dependence, and underlying stochastic processes can each be
substituted. Fitting may also be performed with respect to the output of
physics-based simulations, thereby establishing an indirect link to physical
parameters while achieving much higher computational speed. With these
considerations, the model represents a flexible foundation for implementing
large-scale neuromorphic simulations that incorporate realistic device
behavior.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

TH wrote the manuscript, performed data analysis, and implemented the model.
AE carried out the ($I,U$) measurement. JN and GM fabricated the ReRAM
devices. RW and DW co-advised the project, and DB conceived of the concept and
was in charge of the project.

## Acknowledgments

The authors thank Thomas Pössinger of RWTH Aachen for illustrating Figs 1 and
2.

## Code Availability

A Julia implementation of the model is available on GitHub
(https://github.com/thennen/StochasticSynapses.jl) and archived in Zenodo
(https://doi.org/10.5281/zenodo.6535411).

## References

  * Burr et al. [2017] Burr GW, Shelby RM, Sebastian A, Kim S, Kim S, Sidler S, et al.  Neuromorphic computing using non-volatile memory.  Adv. Phys. X 2 (2017) 89–124.  10.1080/23746149.2016.1259585.
  * Sangwan and Hersam [2020] Sangwan VK, Hersam MC.  Neuromorphic nanoelectronic materials.  Nat. Nanotechnol. 15 (2020) 517–528.  10.1038/s41565-020-0647-z.
  * Milo et al. [2020] Milo V, Malavena G, Monzio Compagnoni C, Ielmini D.  Memristive and CMOS Devices for Neuromorphic Computing.  Materials 13 (2020) 166.  10.3390/ma13010166.
  * Chen et al. [2014] Chen A, Hutchby J, Zhirnov VV, Bourianoff G, editors.  Emerging Nanoelectronic Devices (Chichester, West Sussex, United Kingdom: John Wiley & Sons Inc) (2014).
  * Yu and Chen [2016] Yu S, Chen PY.  Emerging Memory Technologies: Recent Trends and Prospects.  IEEE Solid-State Circuits Mag. 8 (2016) 43–56.  10.1109/MSSC.2016.2546199.
  * You Zhou and Ramanathan [2015] You Zhou, Ramanathan S.  Mott Memory and Neuromorphic Devices.  Proc. IEEE 103 (2015) 1289–1310.  10.1109/JPROC.2015.2431914.
  * Liu et al. [2014] Liu H, Bedau D, Sun J, Mangin S, Fullerton E, Katine J, et al.  Dynamics of spin torque switching in all-perpendicular spin valve nanopillars.  Journal of Magnetism and Magnetic Materials 358–359 (2014) 233–258.  10.1016/j.jmmm.2014.01.061.
  * Waser et al. [2009] Waser R, Dittmann R, Staikov G, Szot K.  Redox-Based Resistive Switching Memories - Nanoionic Mechanisms, Prospects, and Challenges.  Adv. Mater. 21 (2009) 2632–2663.  10.1002/adma.200900375.
  * Kim et al. [2018] Kim S, Lim M, Kim Y, Kim HD, Choi SJ.  Impact of Synaptic Device Variations on Pattern Recognition Accuracy in a Hardware Neural Network.  Sci Rep 8 (2018) 2638.  10.1038/s41598-018-21057-x.
  * Brown et al. [2020] Brown T, Mann B, Ryder N, Subbiah M, Kaplan JD, Dhariwal P, et al.  Language models are few-shot learners.  Larochelle H, Ranzato M, Hadsell R, Balcan MF, Lin H, editors, Advances in Neural Information Processing Systems (Curran Associates, Inc.) (2020), vol. 33, 1877–1901.
  * Jiang and Stewart [2017] Jiang H, Stewart DA.  Using Dopants to Tune Oxygen Vacancy Formation in Transition Metal Oxide Resistive Memory.  ACS Appl. Mater. Interfaces 9 (2017) 16296–16304.  10.1021/acsami.7b00139.
  * Stewart [2019] Stewart DA.  Diffusion of oxygen in amorphous tantalum oxide.  Phys. Rev. Materials 3 (2019) 055605.  10.1103/PhysRevMaterials.3.055605.
  * Kopperberg et al. [2021] Kopperberg N, Wiefels S, Liberda S, Waser R, Menzel S.  A Consistent Model for Short-Term Instability and Long-Term Retention in Filamentary Oxide-Based Memristive Devices.  ACS Appl. Mater. Interfaces 13 (2021) 58066–58075.  10.1021/acsami.1c14667.
  * Ascoli et al. [2015] Ascoli A, Tetzlaff R, Biolek Z, Kolka Z, Biolkova V, Biolek D.  The Art of Finding Accurate Memristor Model Solutions.  IEEE J. Emerg. Sel. Topics Circuits Syst. 5 (2015) 133–142.  10.1109/JETCAS.2015.2426493.
  * Messaris et al. [2018] Messaris I, Serb A, Stathopoulos S, Khiat A, Nikolaidis S, Prodromakis T.  A Data-Driven Verilog-A ReRAM Model.  IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 37 (2018) 3151–3162.  10.1109/TCAD.2018.2791468.
  * Ielmini and Milo [2017] Ielmini D, Milo V.  Physics-based modeling approaches of resistive switching devices for memory and in-memory computing applications.  J. Comput. Electron. 16 (2017) 1121–1143.  10.1007/s10825-017-1101-9.
  * Abbaspour et al. [2020] Abbaspour E, Menzel S, Jungemann C.  Studying the switching variability in redox-based resistive switching devices.  J. Comput. Electron. 19 (2020) 1426–1432.  10.1007/s10825-020-01537-y.
  * Reuben et al. [2019] Reuben J, Fey D, Wenger C.  A Modeling Methodology for Resistive RAM Based on Stanford-PKU Model With Extended Multilevel Capability.  IEEE Trans. Nanotechnology 18 (2019) 647–656.  10.1109/TNANO.2019.2922838.
  * Mayer et al. [2010] Mayer J, Khairy K, Howard J.  Drawing an elephant with four complex parameters.  American Journal of Physics 78 (2010) 648–649.  10.1119/1.3254017.
  * Li et al. [2017] Li H, Huang P, Gao B, Liu X, Kang J, Philip Wong HS.  Device and Circuit Interaction Analysis of Stochastic Behaviors in Cross-Point RRAM Arrays.  IEEE Trans. Electron Devices 64 (2017) 4928–4936.  10.1109/TED.2017.2766046.
  * Maria Puglisi et al. [2015] Maria Puglisi F, Larcher L, Padovani A, Pavan P.  Bipolar Resistive RAM Based on HfO2: Physics, Compact Modeling, and Variability Control.  IEEE J. Emerg. Sel. Topics Circuits Syst. (2015).
  * Bengel et al. [2020] Bengel C, Siemon A, Cuppers F, Hoffmann-Eifert S, Hardtdegen A, von Witzleben M, et al.  Variability-Aware Modeling of Filamentary Oxide-Based Bipolar Resistive Switching Cells Using SPICE Level Compact Models.  IEEE Trans. Circuits Syst. Regul. Pap. (2020) 1–13.  10.1109/TCSI.2020.3018502.
  * Jiang et al. [2016] Jiang Z, Wu Y, Yu S, Yang L, Song K, Karim Z, et al.  A Compact Model for Metal–Oxide Resistive Random Access Memory With Experiment Verification.  IEEE Trans. Electron Devices 63 (2016) 1884–1892.  10.1109/TED.2016.2545412.
  * Chen and Yu [2015] Chen PY, Yu S.  Compact Modeling of RRAM Devices and Its Applications in 1T1R and 1S1R Array Design.  IEEE Trans. Electron Devices 62 (2015) 4022–4028.  10.1109/TED.2015.2492421.
  * Siemon et al. [2019] Siemon A, Wouters D, Hamdioui S, Menzel S.  Memristive Device Modeling and Circuit Design Exploration for Computation-in-Memory.  2019 IEEE International Symposium on Circuits and Systems (ISCAS) (Sapporo, Japan: IEEE) (2019), 1–5.  10.1109/ISCAS.2019.8702600.
  * Bocquet et al. [2014] Bocquet M, Aziza H, Zhao W, Zhang Y, Onkaraiah S, Muller C, et al.  Compact Modeling Solutions for Oxide-Based Resistive Switching Memories (OxRAM).  JLPEA 4 (2014) 1–14.  10.3390/jlpea4010001.
  * Huang et al. [2017] Huang P, Zhu D, Chen S, Zhou Z, Chen Z, Gao B, et al.  Compact Model of HfOx-Based Electronic Synaptic Devices for Neuromorphic Computing.  IEEE Trans. Electron Devices 64 (2017) 614–621.  10.1109/TED.2016.2643162.
  * Wald and Kvatinsky [2019] Wald N, Kvatinsky S.  Understanding the influence of device, circuit and environmental variations on real processing in memristive memory using Memristor Aided Logic.  Microelectronics Journal 86 (2019) 22–33.  10.1016/j.mejo.2019.02.013.
  * Yang et al. [2008] Yang JJ, Pickett MD, Li X, Ohlberg DAA, Stewart DR, Williams RS.  Memristive switching mechanism for metal/oxide/metal nanodevices.  Nature Nanotech 3 (2008) 429–433.  10.1038/nnano.2008.160.
  * Hennen et al. [2021] Hennen T, Wichmann E, Elias A, Lille J, Mosendz O, Waser R, et al.  Current-limiting amplifier for high speed measurement of resistive switching data.  Rev. Sci. Instrum. 92 (2021) 054701.
  * Nail et al. [2016] Nail C, Molas G, Blaise P, Piccolboni G, Sklenard B, Cagli C, et al.  Understanding RRAM endurance, retention and window margin trade-off using experimental results and simulations.  2016 IEEE International Electron Devices Meeting (IEDM) (San Francisco, CA, USA: IEEE) (2016), 4.5.1–4.5.4.  10.1109/IEDM.2016.7838346.
  * Cario and Nelson [1996] Cario MC, Nelson BL.  Autoregressive to anything: Time-series input processes for simulation.  Operations Research Letters 19 (1996) 51–58.  10.1016/0167-6377(96)00017-X.
  * Rezende and Mohamed [2015] Rezende DJ, Mohamed S.  Variational Inference with Normalizing Flows.  Proceedings of the 32nd International Conference on Machine Learning (PMLR) (2015), vol. 37, 1530–1538.
  * Chen [2015] Chen A.  Utilizing the Variability of Resistive Random Access Memory to Implement Reconfigurable Physical Unclonable Functions.  IEEE Electron Device Lett. 36 (2015) 138–140.  10.1109/LED.2014.2385870.
  * Hamilton [1994] Hamilton JD.  Time Series Analysis (Princeton, N.J: Princeton University Press) (1994).
  * Lütkepohl [2005] Lütkepohl H.  New Introduction to Multiple Time Series Analysis (Berlin: New York : Springer) (2005).
  * Fantini et al. [2015] Fantini A, Gorine G, Degraeve R, Goux L, Chen CY, Redolfi A, et al.  Intrinsic program instability in HfO2 RRAM and consequences on program algorithms.  2015 IEEE International Electron Devices Meeting (IEDM) (Washington, DC, USA: IEEE) (2015), 7.5.1–7.5.4.  10.1109/IEDM.2015.7409648.
  * Roldán et al. [2019] Roldán JB, Alonso FJ, Aguilera AM, Maldonado D, Lanza M.  Time series statistical analysis: A powerful tool to evaluate the variability of resistive switching memories.  Journal of Applied Physics 125 (2019) 174504.  10.1063/1.5079409.
  * Seabold and Perktold [2010] Seabold S, Perktold J.  Statsmodels: Econometric and Statistical Modeling with Python.  Python in Science Conference (Austin, Texas) (2010), 92–96.  10.25080/Majora-92bf1922-011.
  * Fantini et al. [2013] Fantini A, Goux L, Degraeve R, Wouters D, Raghavan N, Kar G, et al.  Intrinsic switching variability in HfO2 RRAM.  2013 5th IEEE International Memory Workshop (Monterey, CA, USA: IEEE) (2013), 30–33.  10.1109/IMW.2013.6582090.
  * Dalgaty et al. [2021] Dalgaty T, Castellani N, Turck C, Harabi KE, Querlioz D, Vianello E.  In situ learning using intrinsic memristor variability via Markov chain Monte Carlo sampling.  Nat Electron 4 (2021) 151–161.  10.1038/s41928-020-00523-3.
  * Butcher et al. [2012] Butcher B, Bersuker G, Young-Fisher KG, Gilmer DC, Kalantarian A, Nishi Y, et al.  Hot Forming to Improve Memory Window and Uniformity of Low-Power HfOx-Based RRAMs.  2012 4th IEEE International Memory Workshop (Milan: IEEE) (2012), 1–4.  10.1109/IMW.2012.6213647.
  * Zhao et al. [2014] Zhao L, Chen HY, Wu SC, Jiang Z, Yu S, Hou TH, et al.  Multi-level control of conductive nano-filament evolution in HfO2 ReRAM by pulse-train operations.  Nanoscale 6 (2014) 5698–5702.  10.1039/C4NR00500G.
  * Moon et al. [2019] Moon J, Ma W, Shin JH, Cai F, Du C, Lee SH, et al.  Temporal data classification and forecasting using a memristor-based reservoir computing system.  Nat. Electron. 2 (2019) 480–487.  10.1038/s41928-019-0313-3.
  * Menzel et al. [2015] Menzel S, Böttger U, Wimmer M, Salinga M.  Physics of the Switching Kinetics in Resistive Memories.  Adv. Funct. Mater. 25 (2015) 6306–6325.  10.1002/adfm.201500825.
  * Ma et al. [2019] Ma W, Chiu PF, Choi WH, Qin M, Bedau D, Lueker-Boden M.  Non-Volatile Memory Array Based Quantization- and Noise-Resilient LSTM Neural Networks.  2019 IEEE International Conference on Rebooting Computing (ICRC) (San Mateo, CA, USA: IEEE) (2019), 1–9.  10.1109/ICRC.2019.8914713.
  * Wiefels et al. [2020] Wiefels S, Bengel C, Kopperberg N, Zhang K, Waser R, Menzel S.  HRS Instability in Oxide-Based Bipolar Resistive Switching Cells.  IEEE Trans. Electron Devices 67 (2020) 4208–4215.  10.1109/TED.2020.3018096.
  * Besard et al. [2019] Besard T, Foket C, De Sutter B.  Effective Extensible Programming: Unleashing Julia on GPUs.  IEEE Trans. Parallel Distrib. Syst. 30 (2019) 827–841.  10.1109/TPDS.2018.2872064.
  * Hennen [2022] Hennen T.  StochasticSynapses.jl.  Zenodo (2022).  10.5281/zenodo.6535411.
  * Pedroni et al. [2019] Pedroni BU, Joshi S, Deiss SR, Sheik S, Detorakis G, Paul S, et al.  Memory-Efficient Synaptic Connectivity for Spike-Timing- Dependent Plasticity.  Front. Neurosci. 13 (2019) 357.  10.3389/fnins.2019.00357.
  * Pedroni et al. [2020] Pedroni BU, Deiss SR, Mysore N, Cauwenberghs G.  Design Principles of Large-Scale Neuromorphic Systems Centered on High Bandwidth Memory.  2020 International Conference on Rebooting Computing (ICRC) (Atlanta, GA, USA: IEEE) (2020), 90–94.  10.1109/ICRC2020.2020.00013.
  * Kantorovich [1960] Kantorovich LV.  Mathematical Methods of Organizing and Planning Production.  Management Science 6 (1960) 366–422.  10.1287/mnsc.6.4.366.
  * Park et al. [2013] Park S, Noh J, Choo Ml, Sheri AM, Chang M, Kim YB, et al.  Nanoscale RRAM-based synaptic electronics: Toward a neuromorphic computing device.  Nanotechnology 24 (2013) 384009.  10.1088/0957-4484/24/38/384009.
  * Ambrogio et al. [2016] Ambrogio S, Balatti S, Milo V, Carboni R, Wang ZQ, Calderoni A, et al.  Neuromorphic Learning and Recognition With One-Transistor-One-Resistor Synapses and Bistable Metal Oxide RRAM.  IEEE Trans. Electron Devices 63 (2016) 1508–1515.  10.1109/TED.2016.2526647.
  * Ielmini [2011] Ielmini D.  Modeling the Universal Set/Reset Characteristics of Bipolar RRAM by Field- and Temperature-Driven Filament Growth.  IEEE Trans. Electron Devices 58 (2011) 4309–4317.  10.1109/TED.2011.2167513.
  * Nardi et al. [2011] Nardi F, Ielmini D, Cagli C, Spiga S, Fanciulli M, Goux L, et al.  Control of filament size and reduction of reset current below 10$\mu$A in NiO resistance switching memories.  Solid-State Electron. 58 (2011) 42–47.  10.1016/j.sse.2010.11.031.
  * Nardi et al. [2012] Nardi F, Larentis S, Balatti S, Gilmer DC, Ielmini D.  Resistive Switching by Voltage-Driven Ion Migration in Bipolar RRAM—Part I: Experimental Study.  IEEE Trans. Electron Devices 59 (2012) 2461–2467.  10.1109/TED.2012.2202319.
  * Nishi et al. [2015] Nishi Y, Fleck K, Böttger U, Waser R, Menzel S.  Effect of RESET Voltage on Distribution of SET Switching Time of Bipolar Resistive Switching in a Tantalum Oxide Thin Film.  IEEE Trans. Electron Devices 62 (2015) 1561–1567.  10.1109/TED.2015.2411748.
  * Kim et al. [2016a] Kim W, Menzel S, Wouters DJ, Guo Y, Robertson J, Roesgen B, et al.  Impact of oxygen exchange reaction at the ohmic interface in Ta2O5-based ReRAM devices.  Nanoscale 8 (2016a) 17774–17781.  10.1039/C6NR03810G.
  * La Torre et al. [2016] La Torre C, Fleck K, Starschich S, Linn E, Waser R, Menzel S.  Dependence of the SET switching variability on the initial state in HfO _x_ -based ReRAM.  Phys. Status Solidi A 213 (2016) 316–319.  10.1002/pssa.201532375.
  * Kim et al. [2016b] Kim KM, Yang JJ, Strachan JP, Grafals EM, Ge N, Melendez ND, et al.  Voltage divider effect for the improvement of variability and endurance of TaOx memristor.  Sci Rep 6 (2016b) 20085.  10.1038/srep20085.
  * Piccolboni et al. [2015] Piccolboni G, Molas G, Portal JM, Coquand R, Bocquet M, Garbin D, et al.  Investigation of the potentialities of Vertical Resistive RAM (VRRAM) for neuromorphic applications.  2015 IEEE International Electron Devices Meeting (IEDM) (Washington, DC, USA: IEEE) (2015), 17.2.1–17.2.4.  10.1109/IEDM.2015.7409717.
  * Ielmini et al. [2011] Ielmini D, Nardi F, Cagli C.  Universal Reset Characteristics of Unipolar and Bipolar Metal-Oxide RRAM.  IEEE Trans. Electron Devices 58 (2011) 3246–3253.  10.1109/TED.2011.2161088.