ReARTSim: an ReRAM ARray Transient Simulator with GPU optimized runtime acceleration

The computational models of ReRAM have evolved along with the advancement of ReRAM devices. A shared feature among these models is the utilization of a state variable, denoted as w(t), to represent the resistance of the device at a given time t. This state variable, with a range between 0.0 and 1.0, serves as a proxy for the doping layer thickness within the ReRAM cell. This thickness dynamically changes in response to applied voltage or current. A higher value of w(t) represents a greater doping thickness, indicating a lower resistance, whereas a smaller w(t) value corresponds to higher resistance [27]. The static I–V relation involving the applied voltage, ReRAM cell current, and the state variable can either be linear:

$$\begin{equation}{R_{{\text{MEM}}}} = { }{R_{{\text{on}}}}\cdot w\left( t \right) + {R_{{\text{off}}}}\cdot\left( {1 - w\left( t \right)} \right)\end{equation} \tag{ 2.1a }$$

$$\begin{equation}{i_{\text{M}}} = \frac{{{V_{\text{M}}}}}{{{R_{{\text{MEM}}}}}}\end{equation} \tag{ 2.1b }$$

in which R_MEM is an equivalent ReRAM resistance, R_on is the on-state resistance (minimum) and R_off denotes the off-state resistance (maximum). Alternatively, this relation could be non-linear, as seen in the model described in [28].

Throughout the evolution of ReRAM modeling, two sets of models have been proposed to describe the behavior of the state variable w under applied voltage and current. The first set of models characterizes the rate of change for w in proportional to cell current:

$$\begin{equation}\frac{{dw}}{{dt}} = \alpha \cdot f\left( w \right)\cdot i\end{equation} \tag{ 2.2 }$$

in which α is a scaling coefficient, f(w) denotes a window function, and i is the ReRAM cell current. The window function f(w) sets rate of state variable change to zero at boundary conditions, when w approaches 0 or 1. Another set of models incorporates a programming threshold voltage, observed in ReRAM switching behavior as reported in [29]. This set of model uses the equation as:

$$\begin{equation}\frac{{dw}}{{dt}} = \,\alpha \cdot f\left( w \right)\cdot g\left( V \right)\end{equation} \tag{ 2.3 }$$

in which α and f(w) retain their meanings from equation (2.2), while g(V) is a nonlinear function of applied voltage V, capable of simulating the threshold voltage effect. Various g(V) models have been proposed to emulate the I–V curves of ReRAM devices based on experimental data.

In this study, we have chosen four example models from literatures (referred to as Cell A to Cell D) that fall within the above two sets. These models have been implemented within ReARTSim framework's device library, and detailed information is provided in table 1. For the models exhibiting a linear I–V relationship, we have selected a fixed R_on value of 100 Ω and R_off value of 10 000 Ω (resulting in R_off /R_on = 100, which is a reasonable ratio) to simplify the MVM weight mapping. The remaining device parameters have been chosen based on values reported in relevant literature or assigned approximate values for the purpose of this study. Note that for Cell D, since w represents undoped layer thickness in the original literature [32], we have made slight modifications to both the original equations and signs of parameters, to ensure consistency in the definition of w with respect to this work. It is important to note that these models have been utilized solely for demonstrating the system's ReRAM crossbar transient simulation capabilities. Users have the flexibility to effortlessly modify these model parameters or incorporate additional device models that align more closely with their specific requirements into the framework.

Table 1. The simulated ReRAM model in this work.

Model	Model equation	I–V relation	Model parameters
Cell A [30]	$\frac{{dw}}{{dt}} = { }\alpha \cdot f\left( w \right) \cdot { }{i_{\text{M}}}$	${i_{\text{M}}} = \frac{{{V_{\text{M}}}}}{{{R_{{\text{MEM}}}}}}$	α = 10 000
			p = 10
	$f\left( w \right) = 1 - {\left( {2w - 1}\right)^{2p}}$
Cell B [31]	$\frac{{dw}}{{dt}} = { }\alpha \cdot f\left( w \right) \cdot { }{i_{\text{M}}}$	${i_{\text{M}}} = \frac{{{V_{\text{M}}}}}{{{R_{{\text{MEM}}}}}}$	α = 10 000
			p = 10
	$f\left( w \right) = 1 - {\left( {w - {\text{step}}\left( { - {i_{\text{M}}}} \right)} \right)^{2p}}$
Cell C [28]	$\frac{{dw}}{{dt}} = \,f\left( w \right) \cdot \,g\left( V \right)$	${i_{\text{M}}} = \left\{ {\begin{array}{*{20}{c}} {{a_1}\cdot w\cdot\text{sinh} \left( {b\cdot V} \right),{ }V \unicode{x2A7E} 0} \\ {{a_2}\cdot w\cdot\text{sinh} \left( {b\cdot V} \right),V \lt 0} \end{array}} \right.$	Vp = 1.5
	$g\left( V \right) = \left\{ {\begin{array}{*{20}{c}} {{A_p}\left( {{e^V} - {e^{{V_p}}}} \right),\,V \gt {V_p}} \\ { - {A_n}\left( {{e^{ - V}} - {e^{{V_n}}}} \right),\,V \lt - {V_n}\,} \\ {0,\, - {V_n} \unicode{x2A7D} V \unicode{x2A7D} {V_p}} \end{array}} \right.$		Vn = 0.5
			Ap = 0.005
			An = 0.08
			a1 = 3.7(10−7)
			a2 = 4.4(10−7)
Cell D [32]	$\frac{{dw}}{{dt}} = { }f\left( w \right) \cdot { }g\left( V \right)$	${i_{\text{M}}} = \frac{{{V_{\text{M}}}}}{{{R_{{\text{MEM}}}}}}$	Von = 0.2
	$\begin{array}{c}\\\\[40pt]g\left( V \right) =\left \{ \begin{gathered} {k_{{\text{off}}}}\left( {{{\left( {\frac{V}{{{V_{{\text{off}}}}}} - 1} \right)}^{{\alpha _{{\text{off}}}}}}} \right),V \lt {V_{{\text{off}}}} \hfill \\ {k_{{\text{on}}}}\left( {{{\left( {\frac{V}{{{V_{{\text{on}}}}}} - 1} \right)}^{{\alpha _{{\text{on}}}}}}} \right),V \gt {V_{{\text{on}}}} \hfill \\ 0,{V_{{\text{off}}}} \unicode{x2A7D} V \unicode{x2A7D} {V_{{\text{on}}}} \hfill \\ \end{gathered} \right.\end{array}$		Voff = −0.02
			kon = 10
			koff = −5(10−4)
			αon = 3

We simulate the time-stamped transient I–V response of the models using the Euler method:

$$\begin{equation}\frac{{dw}}{{dt}}\left( {{t_0}} \right) \approx \frac{{w\left( {{t_0} + h} \right) - w\left( {{t_0}} \right)}}{h}\end{equation} \tag{ 2.4 }$$

a fundamental numerical method to solve differential equations. In this method, the simulation timestep size, denoted as h, can have the unit of nanoseconds (ns), microseconds (µs), or milliseconds (ms). The variable w(t₀) denotes the known state variable at time t₀ , and w(t₀ + h) represents the state variable at the subsequent timestamp t₀ + h. By utilizing the calculated $\frac{{dw}}{{dt}}\left( {{t_0}} \right)$ obtained from equation (2.2) or equation (2.3) and substituting it into equation (2.4), the time-stamped function w(t) and consequently the transient I–V curve can be derived. More advanced numerical methods for solving differential equations, such as the Runge–Kutta method, can be implemented in the future to enhance simulation accuracy [33]. In the example simulations, we applied triangular input voltage waves with both positive and negative amplitudes. Figures 2(a), (c) and (e) depict the simulated time-stamped voltage and current responses of Cell B, C and D. Figures 2(b), (d) and (f), offering an alternative perspective on figures 2(a), (c) and (e), illustrate the transition in cell current under the applied voltage V, with arrowheads indicating the rising and falling slopes of the applied triangular voltage. The triangular pulse counts are labeled for these plots. The simulated curves align with those reported in the initial paper [28, 31, 32]. As further detailed in section 2.2, we leveraged the Euler method along with the device model outlined in table 1 to conduct comprehensive transient simulations of the ReRAM crossbar.

Zoom In Zoom Out Reset image size

As introduced in the earlier sections, an array of ReRAM cells can be utilized to construct a ReRAM crossbar, and used for performing MVM computation. Figure 3 provides an illustration of a typical ReRAM crossbar structure. Each individual ReRAM device possesses a conductance value denoted as G_rc and is linked to a wordline in its corresponding row and a bitline in its column. Upon applying a voltage V_r to the wordline in row r, the current in the bitline at column c will either be V_r × G_rc if the column is activated, or 0 if the column is deactivated. Due to finite values of R_on and R_off, the conductance G_rc is confined within the range [1/R_off, 1/R_on]. In our ReRAM model, where we have chosen fixed R_on and R_off values (100 Ω and 10 000 Ω), the achievable conductance range within the crossbar spans from 0.0001 to 0.01 S. For standardization, we have chosen a conductance value of 0.001 S as our base unit, with this value corresponding to a matrix value of 1.0. Consequently, valid matrix values achievable by the ReRAM crossbar span from 0.1 to 10, covering common matrix magnitude. Given that the current is summed along the column, in the case of a typical right MVM, the original matrix values along each row correspond to the ReRAM crossbar cell conductance along each column. As an illustration, consider a 2 × 2 matrix with values of [[0.2, 0.3], [0.8, 0.6]]; it will be mapped onto a 2 × 2 ReRAM crossbar with conductance values of [[0.0002, 0.0003], [0.0008, 0.0006]], where [0.0002, 0.0003] is mapped along the first column, and [0.0008, 0.0006] along the second column. These values fall within the achievable range of ReRAM cell conductance values.

Zoom In Zoom Out Reset image size

In the process of MVM computation, the input vector is converted from digital to analog voltages in the unit of volt (V) using a digital-to-analog converter (DAC), subsequently driving the respective wordlines. For instance, an input vector of [1, 1.5] will be transformed into a 1 V pulse on the first row and a 1.5 V pulse on the second row, each with a specific duration. During the MVM computation, all columns are activated simultaneously, and within each bitline, the currents from all ReRAM devices in the same column are summed. The bitline currents are sampled by a sense amplifier (SA) and analog-to-digital converter, transforming them into digital outputs that constitute the output vector. As the base unit of conductance is set to 0.001 S, the bitline current value, measured in mA, is directly utilized as the output vector element. Continuing from the previous example, the output current in the first column is computed as 1 × 0.0002 + 1.5 × 0.0003 = 0.00065 A, which translates to 0.65 mA. This value corresponds to the first element of the MVM output vector, aligning with the direct MVM computation result. Notably, parasitic trace resistance and capacitance are present within the crossbar, exerting an impact on the computation outcome, particularly at high speeds. Streaming data can be fed into the input processing logic and results are retrieved at the output data processing logic, which are coordinated and controlled by the central control logic.

Several non-idealities could exist in ReRAM crossbar arrays. It is established that ReRAM switching behavior is subject to variations from device-to-device and cycle-to-cycle [34, 35]. Device-to-device variation is typically influenced by non-ideal fabrication processes [36], while cycle-to-cycle variation may arise due to stochastic effects during filament formation in ReRAM resistive switching [37]. Such variations can potentially lead to incorrect computation results in MVM. Proposed training methods, such as those outlined in [38, 39], seek to mitigate these issues by considering device variations during ReRAM crossbar network training. Meanwhile, similar to any resistive component, thermal noise can contribute to output bitline current. As a result, it is important for the software to possess the capability to conduct simulations while accounting for these non-idealities.

The ReARTSim simulation framework is implemented using the C++ programming language, featuring a class-based structure depicted in figure 4(a). The Cell class serves to describe the behavioral attributes of the ReRAM cell. Within this class, the ReRAM cell's model, model parameters, state variable w, and time-stamped voltage, current, and power are included. It also includes values for parasitic trace resistance and capacitance, essential for high-speed simulation. The initiation of time-stamped voltage originates from user-defined inputs, and one approach is through the use of the Instruction class, detailed below. Subsequently, by executing the runSimEuler() method, the temporal evolution of w is computed in accordance with equations (2.2)–(2.4), simultaneously updating the values of current and power.

Zoom In Zoom Out Reset image size

The Instruction class is utilized to generate a set of user-specific instructions, providing guidance on supplying time-stamped voltage input to individual ReRAM cells. These instructions cover a variety of operational scenarios, spanning basic cell read/write operations as well as more intricate MVM instructions. In the case of the read/write instruction, which targets a single cell, voltage is applied to a specific wordline for a specified duration, and the resulting current is collected from the corresponding activated bitline. Conversely, the MVM instruction involves applying voltage across all wordlines based on the input vector for a specified duration, followed by the collection of currents from all bitlines. The duration of applied voltage in read/write or MVM instruction is a variable specified by user.

The Array class functions as the foundational component that defines a ReRAM crossbar array. It incorporates a 2D vector containing Cell objects that represent ReRAM cells within the crossbar array. Additionally, it includes an array of Instruction objects, facilitating user input. Time-stamped data, including voltage for each row (VvDrives), the active status for each column (VvActives), current across columns (IvSenses), and the aggregate crossbar power (PvSys) are also stored within this class. Simulation parameters, including the timestamp unit (ns, µs, or ms), as well as the total count of timestamps in the simulation, can be configured in the Array class. For simulations requiring an extensive timestamp count, the software optimizes CPU memory usage by segmenting the simulation into multiple cycles. Within its runSimEuler() method (equally named at both array and cell levels), the Array object initiates by processing the instruction list in the current simulation cycle, subsequently updating the VvDrives and VvActives vectors. Following this, each individual cell invokes its respective runSimEuler() methods to simulate cell responses within the CPU-based workflow. Upon completion of cell simulations, the Array class computes the IvSenses and PvSys vectors by summing values from individual cells.

The Matrix class is developed to facilitate the initialization of ReRAM cells for subsequent use in MVM operations. This class contains a 2D vector that stores the original weight matrix values, which indicate the desired conductance for each ReRAM cell in the crossbar. In the Matrix class, a MatrixHelper object is employed to assist in generating write instructions for initializing the conductance of ReRAM cells. Within the MatrixHelper object, a pre-calculated data structure, gmap, is constructed in advance for each ReRAM cell model. This gmap maps a finite number of ReRAM conductance values to specific write instructions with fixed applied voltage levels and durations. To set the ReRAM cell to its target conductance value, the Matrix object queries the MatrixHelper to search the pre-calculated gmap, determining the optimal write instruction. In our framework, we have implemented a basic weight initialization algorithm wherein cells are individually programmed, which is briefly discussed in section 3.1. In this manner, for an r × c ReRAM array, the Matrix object generates a total of r × c write instructions. More efficient crossbar initialization algorithms, as described in [40], can also be incorporated into the framework.

To facilitate the need for verifying MVM results, the Array class includes functions capable of sampling the bitline current for each MVM operation. Subsequently, these functions generate an output vector using the bitline current values measured in mA, which are then compared with reference values computed directly through mathematical equation. The sampling takes place at a user-defined delay after the initiation of each MVM instruction. Additionally, a ConvKernel class has been developed to support 2D matrix convolution operations on the ReRAM crossbar, detailed further in section 3.2. In brief, the 2D convolution is translated into a sequence of MVM operations, utilizing a fixed weight matrix that has been programmed in the crossbar. The Array object then samples the bitline current and constructs convolutional output for this sequence of MVM operations, facilitating a comparison with the computed reference.

While the primary focus of the framework is on cell-level transient simulation, it also offers performance and power benchmarks, similar to non-SPICE-based ReRAM simulators. Following the procedures outlined in [22], the array object can report performance in giga floating-point operations per second (GFLOPS), energy and power in nJ and W, and energy efficiency in GFLOPS/W after completing each batch of MVM simulations. Users can leverage these metrics to analyze system power efficiency and determine the peak MVM performance for each array and cell model before any inference accuracy drop occurs. A more detailed example is provided in section 3.3, the AI inference section.

To account for device-to-device variations, the software possesses the capability to import pre-defined device parameters into each ReRAM Cell object. Furthermore, APIs have been provided to enable users to configure ReRAM cell parameters with specific trends of change, such as linear changes along rows or columns, aiming to mirror actual non-idealities encountered during the fabrication process. The software also permits users to set the noise level across bitlines within the simulation. These features facilitate the verification of training methods proposed in [38, 39] for accommodating device variations, and [41] for enhancing noise tolerance, utilizing the framework.

The flowchart for a typical MVM computation program performed in this framework is depicted in figure 4(b). Upon simulation initiation, users begin by utilizing APIs to import both the weight matrix and a list of input vectors. A pre-processing step is executed on the imported weight matrix to constrain the matrix values within the range of 0.1 and 10, suitable for mapping onto the ReRAM crossbar. To initialize ReRAM cells, a Matrix object based on the imported weight matrix is declared first. Then APIs from the Matrix class are utilized to generate a specific write instruction for each weight matrix element, facilitating the setting of cell conductance to the desired value. Additionally, APIs have been developed to generate an MVM instruction for each input vector, where the voltage on each wordline aligns with the corresponding input vector values. The duration of the MVM instruction can also be specified through the APIs. Users then instantiate an Array object, specifying the ReRAM cell type and optionally incorporating ReRAM cell parameters. Throughout the simulation phase, users can select from two workflow choices: the CPU-based workflow and the GPU-based workflow. The CPU-based workflow, devoid of the need for specialized hardware support, conducts transient simulations one cell after another, in line with previous descriptions. The GPU-based workflow, that will be described in detail in section 2.3, requires GPU hardware installed to accelerate computation. In the GPU-based workflow, the CPU loads the simulation data—including cell parameters and timestamped row voltages—onto the GPU global memory, then waits for GPU-driven crossbar simulations. Subsequent to simulation, the CPU retrieves the data sent back from GPU. The CPU then undertakes data processing, perform tasks including bitline sampling to derive MVM outputs, and subsequent comparison against reference values.

The GPU, initially designed for display image processing, has gained attraction in various computational fields due to its architecture featuring numerous parallel processing cores [42]. In the context of large ReRAM crossbar arrays, the computation of transient responses for a total of r × c cells (where r represents row count and c denotes column count) through a single-threaded program can prove time-consuming. While multithreaded CPU programming might enhance performance, scalability is restricted by the number of available CPU cores. To address this limitation, we have developed GPU kernel software utilizing the Nvidia CUDA library [26] to accelerate ReRAM simulations. Utilizing the common structure of each ReRAM cell in the crossbar, the GPU emerges as a robust hardware choice for simulating crossbar arrays. Its capability to handle programming patterns involving thousands of threads performing similar tasks reinforces its suitability for this purpose.

Figure 5 illustrates the architecture of the GPU acceleration software deployed in the GPU-based workflow. Analogous to the CPU-based workflow, the entire simulation duration is partitioned into multiple simulation cycles. Within each cycle, the CPU initiates the process by requesting the GPU to allocate global memory for all cells, which is used to store cell parameters, state variable w, time-stamped voltage, and current. These values are then loaded onto the GPU's global memory. Subsequently, a total of r × c computation kernels are instantiated on the GPU. These kernels, akin to the runSimEuler() methods in CPU-based workflow, calculate the transient response for each of the r × c cells. After these computation kernels are completed and synchronized, reduction kernels are used to aggregate current across each bitline and overall array power. With the simulation time in one cycle partitioned into t segments, a total of t reduction kernels are launched, with each reduction kernel computing current and power for one such segment. After the completion and synchronization of the reduction kernels, the GPU forwards the computed transient data to the CPU for post-processing, awaiting the next simulation cycle computation. To further reduce simulation time, user can opt to stop the retrieval of individual cell simulation result from the GPU. When this option is activated, only bitline current and total power are transferred back to CPU memory post-simulation. which cuts down the data transfer time. Given that modern GPUs can accommodate thousands of threads and blocks, substantial parallelization of these kernels can be achieved [43, 44].

Zoom In Zoom Out Reset image size

To assess the transient simulation capabilities of the ReARTSim framework, we executed illustrative tasks including matrix weight initialization and MVM computations. The operation schematics for these tasks are depicted in figures 6(a) and (b), respectively. In the matrix weight initialization task, a 4 × 4 ReRAM Array object with ReRAM cell type C was declared, chosen for its simplicity and inclusion of threshold programming voltage. It was slightly adjusted to incorporate a linear I–V relation, resulting in a fixed resistance under varying applied voltages. This adjustment was made to demonstrate the framework's MVM capabilities. Additionally, a 4 × 4 Matrix object was declared, representing the target weight matrix to be programmed onto the ReRAM array. The matrix element values were constrained within positive ranges, allowing their conversion to valid ReRAM cell conductance. Following this, the Matrix object generated a series of write instructions to initialize ReRAM cells. As described earlier, we implemented a basic weight initialization algorithm, programming each cell individually. The operation is illustrated in figure 6(a), where the voltage on each wordline was applied until all cells in the present row were initialized, and then the voltage was shifted to the next wordline. Bitlines were activated in a round-robin fashion, allowing cells from different columns to be programmed consecutively. The duration that cells in row m were being programmed was labeled as r_m , and this duration is marked in figure 6(a).

Zoom In Zoom Out Reset image size

The simulation was then executed until all ReRAM cells within the crossbar were initialized to the desired target weight. Figure 7(a) presents the transient current across the four bitlines, while figure 7(b) displays the transient power of the array. As seen from the figures, the time scale for the cell weight initialization task is relatively lengthy (measured in seconds). This extended duration results from the time required to alter the cell's state variable w, which is in line with the ReRAM cell model. The figures show noticeable peaks in current and power (∼1.5 s). Additionally, the user API supports the printing of the initialized conductance in comparison to the reference matrix. This functionality facilitates the evaluation of the accuracy of the matrix weight initialization process.

Zoom In Zoom Out Reset image size

In the MVM computation task, we imported multiple input vectors of length 4, declaring MVM instructions for these vectors, and performed transient simulations on the Array object. As illustrated in figure 6(b), during each MVM operation, the voltages converted from the input vector were simultaneously applied to all wordlines, and all bitlines were activated concurrently. The current from each bitline was then amplified by SA and then subsequently sampled by DAC to form the output vector. Because MVM computation can be significantly faster compared with weight initialization task, since the state variable w change is not involved, a shorter simulation timescale (ns) was used for plotting the transient response. Additionally, parasitic resistors and capacitors were enabled to enhance the precision of modeling. Figure 7(c) shows the transient current in the four bitlines. The figure shows the effect of the first-order RC response due to the trace resistance and capacitance, with sampling conducted when all bitline currents are stable. Figure 7(d) presents the console output of the simulated array's MVM result, along with the reference MVM result computed directly. The minor discrepancy between the simulated MVM and the golden result is due to the deviations during the initialization of the weight matrix. Other factors, including device variations, voltage and current noise, and sampling time, may also contribute to minor computational discrepancies and can be assessed through the console output.

To assess the performance enhancement achieved through the transition from CPU-based to GPU-based workflows, we conducted transient simulations for MVM tasks across various ReRAM crossbar array sizes using both CPU and GPU platform. We measured the computation time for each task, and compared the values between two workflows. Figure 8 shows the MVM computation time for array sizes of 8 × 8, 16 × 16, 32 × 32, 64 × 64, and 128 × 128, all including a total of 10 000 timestamps, for both CPU and GPU platform (64 × 64 and 128 × 128 datapoints are indicated in the text rather than plotted). The findings reveal that GPU computation outperforms CPU for larger crossbar sizes—specifically, arrays larger than 16 × 16. At array sizes of 64 × 64 and 128 × 128, GPU performance is approximately 50 times faster than CPU. This disparity can be explained that when array size is small, the CPU is faster due to absence of the GPU kernel launch overhead, and relatively limited parallelism is achievable by GPU kernels. As the number of cells running in parallel increases, the kernel launch overhead becomes insignificant, and GPU kernels leverage higher parallelism [42, 45]. Still, the performance is constrained to a certain limit, which is likely due to GPU memory bound since only global memory is employed in the GPU acceleration kernel [46]. Further optimization strategies, such as utilizing shared memory and wrapped synchronization, could be adopted to refine performance [42, 47].

Zoom In Zoom Out Reset image size

To exemplify the simulation capabilities of the ReARTSim framework for practical applications on the ReRAM crossbar platform, we selected sample NN architectures. We integrated high-level APIs to handle the neural network structures and inputs, and deployed part of their computations onto our simulation framework. Figure 9(a) shows two representative NN architectures: LeNet [48] and Alexnet [49], that were simulated in our framework. Both architectures feature common computation layers: convolution layer (CONV) and fully-connected layer (FC). The FC layer includes MVM with non-linear activation for the output vector, and the MVM is a direct match to the MVM instruction in our framework. However, post-training, matrix weight values typically span from −1.0 to 1.0, while our ReRAM crossbar valid weight lies between 0.1 and 10. As a result, for MVM in the FC layer, we employed a differential-matrix method to handle weight values beyond the intended range. Two auxiliary matrices—each containing positive and negative values extracted from the original matrix—were independently constructed. A uniform offset was then applied to both matrices, elevating all weight values above 0.1 to address elements within the 0–0.1 absolute value range. After MVM computation, the output vector from the negative matrix was subtracted from the positive matrix output vector to yield the final MVM output.

Zoom In Zoom Out Reset image size

For the CONV layer, we employed the technique outlined in [50] to transform the convolution computation into a series of MVM operations. The conceptual process is depicted in figure 9(b). The computation kernel, characterized by a 4D dimension of C_out × C_in × h × w (where C_out is the number of output channels, C_in is the number of input channels, and w and h are the 2D convolution kernel size), was flattened into a 2D vector of dimensions (C_in × h × w) × C_out. Similar to the convolution operation, an input vector of size C_in × h × w was retrieved from the input data at each convolution step, with a sliding window traversing along h and w directions for subsequent convolution steps. This input vector was multiplied with the flattened matrix, equivalent to the element-wise convolution operation. The resulting output vector from the MVM operation was used to update the convolution output at a specific location (h_i, w_j ) across all C_out output channels in the convolution output layer. In our framework, we have introduced a ConvKernel class responsible for holding the convolution kernel. This class features a flatten() method to transform the 4D convolution kernel into a 2D matrix, and can help to compute the convolution reference result. Padding zeros and variable stride lengths are supported in the framework. The Array class incorporates methods to store the sampled MVM outputs from the convolution steps, and updates the outcomes when a convolution operation concludes. Similarly, the previously mentioned differential matrix approach used in MVM was employed for convolutions.

Finally, we integrated our framework with the PyTorch API and conducted complete neural network inference. Initially, we trained the LeNet model using the CIFAR-10 dataset and utilized the built-in pre-trained model for AlexNet. We then input the testing data into both LeNet and AlexNet models and analyzed the inference result. During the inference process, the first two convolutional layers from both models were executed within our ReRAM crossbar simulation framework (highlighted by the red dashed rectangular in figure 9(a)) and the remaining layers were computed using PyTorch. The performance and power benchmarks of the first CONV layer for both NN network architectures, including total simulation time in ns, total number of floating-point (FP) operations, performance measured in GFLOPS, overall energy and power quantified in nJ and W, and energy efficiency measured in GFLOPS/W, are shown in figure 9(c). Similar performance benchmarks for the other NN layers can be presented in a comparable manner. Figure 9(d) shows the inference accuracy under different NN throughput of the first CONV layer in both NN networks, which reflects different MVM computation speed. At higher MVM speed, the first-order RC response is more apparent, and the sampling is more likely to occur before the bitline current stabilizes, leading to a drop in inference accuracy. Figure 9(e) shows the inference accuracy (top-1 error for LeNet, top-5 error for AlexNet) under varying noise levels, while figure 9(f) shows the inference accuracy with different levels of device variation. Since we trained the NN architecture in the conventional approach, a decline in inference accuracy is expected at higher noise and device variation level. This serves to demonstrate the functionality of our simulator.

In addition to its application in neural network inference, our framework extends its utility through the development of APIs capable of simulating other computation tasks that can be effectively mapped onto the ReRAM crossbar network. Notably, we implemented a simulation of principal component analysis (PCA), a statistical technique used for dimensionality reduction in input datasets based on computed feature vectors. Traditionally, PCA involves computationally intensive processes such as computing the covariance matrix and solving for eigenvalues and eigenvectors of that covariance matrix, with costs increasing quadratically with input dimensions [51]. Alternatively, we employed Sanger's Rule, also known as the Generalized Hebbian Algorithm [52], to compute principal components. The algorithm suggests that for an input data dimension of n, we initialize an m × n matrix with random values. Weight updates for each input data point are calculated using the equation:

$$\begin{equation}\Delta {w_{ij}} = \,\alpha \left( {{x_j}{y_i} - {y_i}\mathop \sum \limits_{k = 1}^i {w_{kj}}{y_k}} \right)\end{equation} \tag{ 3.1 }$$

in which α is the learning rate, x_j is the input at index j, y_i is the output at index i, and w_ij represents the weight in row i and column j of the matrix. Initially, vector y is computed by multiplying input data x with the matrix. Subsequent steps involve using y as input and transposing the matrix for further computation $\mathop \sum \nolimits_{k = 1}^i {w_{kj}}{y_k}$. Through iterative processing of input data, the m rows of the initial m × n matrix progressively approach convergence with the first m principal components.

The generalized Hebbian algorithm has been successfully implemented in a ReRAM crossbar, as demonstrated in [53]. Figure 10(a) demonstrates the custom API flow adopted within this software framework. The initial m × n matrix is mapped onto the crossbar with n rows and m columns, following the previously described MVM mapping routine. During the initial timestamp, the input vector x is applied across the rows, and the resulting current across the columns is sampled to compute the output vector y. This y vector is stored, then with all elements marked as zero except for the first element, y[0]. In the subsequent m timestamps, the y vector, which is progressively updated, is applied across the columns, transitioning one value at a time from zero to the previously stored value. The transition occurs incrementally, start with y[1], followed by y[2], and so on, in accordance with the term $\mathop \sum \nolimits_{k = 1}^i {w_{kj}}{y_k}$ defined in equation (3.1). Current sampling takes place across the rows for each timestamp. In the Instruction class, a backward MVM instruction is introduced to facilitate applying the input voltage across the columns and sampling output current at each row within the crossbar. This PCA concept is effectively illustrated using 2D Gaussian data. In figure 10(b), the sampled current across columns and rows during a single computation step is depicted, while figure 10(c) illustrates the computed principal components using the framework.

Zoom In Zoom Out Reset image size