IR-QNN Framework: An IR Drop-Aware Offline Training of Quantized Crossbar Arrays

Resistive Crossbar Arrays present an elegant implementation solution for Deep Neural Networks acceleration. The Matrix-Vector Multiplication, which is the corner-stone of DNNs, is carried out in <inline-formula> <tex-math notation="LaTeX">$O(1)$ </tex-math></inline-formula> compared to <inline-formula> <tex-math notation="LaTeX">$O(N^{2})$ </tex-math></inline-formula> steps for digital realizations of <inline-formula> <tex-math notation="LaTeX">$O(log_{2}(N))$ </tex-math></inline-formula> steps for in-memory associative processors. However, the IR drop problem, caused by the inevitable interconnect wire resistance in RCAs remains a daunting challenge. In this article, we propose a fast and efficient training and validation framework to incorporate the wire resistance in Quantized DNNs, without the need for computationally extensive SPICE simulations during the training process. A fabricated four-bit Au/Al<sub>2</sub>O<sub>3</sub>/HfO<sub>2</sub>/TiN device is modelled and used within the framework with two-mapping schemes to realize the quantized weights. Efficient system-level IR-drop estimation methods are used to accelerate training. SPICE validation results show the effectiveness of the proposed method to capture the IR drop problem achieving the baseline accuracy with a 2% and 4% drop in the worst-case scenario for MNIST dataset on multilayer perceptron network and CIFAR 10 dataset on modified VGG and AlexNet networks, respectively. Other nonidealities, such as stuck-at fault defects, variability, and aging, are studied. Finally, the design considerations of the neuronal and the driver circuits are discussed.


I. INTRODUCTION
Artificial Intelligence hardware acceleration has attracted significant interest [1], [2] especially accelerating deep neural networks (DNNs) with in-memory processing, alleviating the memory-wall bottleneck problem in the Von-Neumann computing architecture. In-memory computing paradigm offers a powerful tool for accelerating artificial intelligence and machine learning algorithms where the matrix-vector multiplication (MVM) computation can be performed in O (1) with resistive crossbar structures and in O(log 2 (N )) with an in-memory associative processor [3], [4], unlike other digital implementations that require O(N 2 ) steps.
Recent advances in non-volatile memory devices, such as Resistive Random Access Memory (RRAM) (memristor), The associate editor coordinating the review of this manuscript and approving it for publication was Shiping Wen .
Spin-transfer torque magnetic random-access memory (STTRAM), and Phase-change memory (PCM) that form the crossbar structure, promise an efficient implementation of MVM computation. The data are stored locally inside the crossbar array eliminating the memory-access need. The ability to efficiently perform MVM computations, especially with RRAM crossbar arrays (RCAs), is very appealing for DNNs since almost all DNN computations can be reduced to MVM operations [5], [6]. Recently, memristor-based crossbar arrays have been used to implement and accelerate many machine learning and deep learning algorithms including pattern recognition [7]- [11], sentimental analysis [12], neural processing [13], reinforcement learning [14] and neuromorphic systems [15].
Recently many RRAM devices have been introduced, experimentally showing the ability to be programmed to multiple distinct states, from 2 states (1-bit) up to 64 states (6-bit) in fabricated devices [16]. These devices have been applied to in-memory computing involving weight parameters for quantized neural networks. Practically, there are multiple technical challenges to realizing RRAM-based DNN hardware. One important issue is how to overcome RRAM arrays' intrinsic inaccuracy and provide an accurate result at the application level, without the need for online (i.e., in-situ) training.RCAs' computation result is susceptible to the device's nonidealities such as device variability, read/write disturbance, stuck-at fault defects, etc. in addition to the inherent noise of analog computation [17], [18]. In particular, the IR drop 1 problem is caused due to the interconnect wire resistance, which affects the MVM computation quite significantly for large arrays and/or advanced technologies [19]- [23]. The interconnect wire resistivity is exponentially increasing with decreasing the technology node due to the increase in electron scattering at grain boundaries, surfaces, and interfaces according to the International Technology Roadmap for Semiconductors (ITRS) roadmap [24]. The intuitive way to mitigate the wire resistance is by improving the metal technology (i.e., less resistivity) or by increasing the interconnect dimensions, which might prevent 3D integration and causes a reduction in density per unit area. For instance, In [25], the authors had to replace nanowire by nano-wall to be able to fabricate a 2nm device to reduce the wire resistance from 10 5 /µm to 10 2 /µm which is not desirable for high dense RCA, especially for 3D integration technologies.
The prior works show that the IR drop problem is a leading cause for degraded performance in MVM computations [26], [27]. For instance, our study based on SPICE-simulationintegrated QNN inference finds that the recognition accuracy for MNIST and CIFAR10 datasets suffers a huge drop as depicted in Fig. 1, which is not acceptable. The higher the wire resistance, the higher accuracy drop.Since SPICE simulation is computationally expensive, it cannot be used for training of QNN. Hardware solutions, such as using 1T1R structure to activate one column at a time, increase the time complexity of MVM to O(N) and require extra hardware to store data [28]. Another solution, such as using 1S1R, helps to mitigate the IR drop problem while achieving the same parallelism of passive structure. However, it highly disturbs the MVM computation due to the exponential nonlinearity of the selector [29].
Prior works typically include numerical or SPICE IR-drop simulations or hardware experiments during training, which highly impacts the training procedure. In [27], the authors proposed an additive noise injection technique, referred to as NIA, to compensate for the effect of the IR drop where the crossbar output current shift with or without IR-drop are collected for all of the testing data. The impact of IR-drop is approximated as a Gaussian noise source at each crossbar output end, with crossbar-wise mean and standard deviation 1 The term IR drop originates from the fact that the amount of voltage drop is proportional to the product of the current (commonly denoted by I) and the wire resistance (indicated by R). extracted from the collected statistics. Finally, the network is retrained with the approximated additive IR drop noise applied at each output of the crossbar array. This technique requires a substantial analog memory (or ADC and memory) to save the analog output currents of each crossbar array, or at least it would involve an iterative SPICE simulation of the entire network, which is not practical for large networks and is not even preferred for small networks. Another method is introduced in [30], where an iterative post-processing technique is proposed for finding the best weight matrix under the IR drop that is very close to the trained weight, which required at least seven iterations for 0.005 MSE. Each iteration requires one IR drop simulation for each weight matrix, which is not practical for large networks. Moreover, If the RCA hardware is available, a software-hardware co-design approach is used where the inference is wholly carried out on the RCA hardware that contains all the nonidealities. The backward path is performed on the software level with ideal parameters [7]. This approach is optimal at the expense of a prolonged training cycle, making it unpractical for DNNs that have billions of parameters, and the quality of model transferability is not guaranteed. In addition, in [31]. we proposed in two neural networks models to model IR drop problem in RCAs. These models were designed and optimized for binary neural networks and require an IR drop dataset to train for each IR drop scenario.
Our work aims to avoid any kind of numerical/SPICE/ Hardware experiments for the IR drop problem during the training and also to avoid any kind of retraining, which would require extra hardware, leading to an increase in the power and area. Our technique directly maps the weights to the hardware without retraining or extra hardware. We use the SPICE-equivalent simulator for validation only, and the proposed techniques are independent of the weight values or network structure. These techniques can also involve hardware for a more accurate estimation of the IR drop problem. The contributions of this work are summarized in the following points: • We first model a fabricated four-bit Au/Al 2 O 3 /HfO 2 /TiN device for accurately mapping the weights to the device's conductances where we introduce two possible mappings for quantized neural network realization.
• We then explain the IR drop problem and evaluate the performance degradation quantitatively.
• An IR-QNN training and validation framework tailored for RRAM crossbar array hardware is introduced, considering the IR drop problem.
• We also introduce efficient software-level methods to incorporate the effect of the IR drop problem without the need for running extensive and time-consuming SPICE simulations.
• We experimentally show that the proposed methods capture the IR drop problem showing a 2% and 4% drop in the worst-case scenario for MNIST dataset on multilayer perceptron network and CIFAR 10 dataset on modified VGG and AlexNet networks, respectively.
• We evaluate the effect of other nonidealities, such as stuck-at fault defects, variability, and retention on the performance.
• Finally, we discuss the design requirements of the neuronal and the driver circuits to guide the designers to have robust and efficient circuits.
This article is organized as follows: Section I discusses the hardware accelerator's hardware realization, where we introduce the multi-bit RRAM device model and the weight mapping. The IR-drop problem is explained in detail in Section III. Section IV discusses the proposed framework for the training and inference and the proposed software methods to estimate the IR drop problem without involving SPICE simulations. Section V discusses the training results and studies the effects of other nonidealities such as device variability, stuck-at fault defects, and aging. Finally, peripheral circuit requirements are discussed in Section VI.

II. IN-MEMORY MVM IMPLEMENTATION
In this section, we discuss the In-Memory MVM using RRAM-crossbar arrays (RCAs), where the weight matrix is stored in the RRAMs. We first model the multi-bit device that is used to realize the weights. Then, we discuss the quantized weight mapping onto the multi-bit device.

A. MULTI-BIT RRAM DEVICE UNDER STUDY
The authors in [16] demonstrated the fabrication of Au/Al 2 O 3 /HfO 2 /TiN RRAM device with a junction area of 10 × 10µm 2 patterned via photolithography and followed by a wet-etching process. The thicknesses of the top electrode (Au) and the bottom electrode (TiN) were 150 and 100 nm, respectively, and the switching material thicknesses are 2 and 6 nm for Al 2 O 3 and HfO 2 , respectively. This device was optimized to have self-compliance and gradual set-switching behavior and is capable of generating up to 16 states with good reliability.
To precisely program the RRAM device, an incremental step pulse programming technique with error correction is used. Starting from a low conductance state, incremental step pulses are applied to the device using a pulse generator until the device reaches the required state. Figures 2a and 2b show the gradual incremental-step programming and the current-voltage hysteresis of the programmed device under different programming conditions. Figures 3a and 3b show the histogram and cumulative distribution function of the measured conductances, respectively. It is worth mentioning that the conductance's variation is due to cycle to cycle read variability. The measured samples are curve-modeled into a Gaussian distribution, and we have found that the mean value of the device's conductance can be modeled as G i = 14+6×i (µS) where G i is the i th highest conductance state for i ∈ [1,15] while the low conductance state is 46.7nS. In order to incorporate the variability in the hardware simulations, there are two ways to integrate the device model in hardware simulation: (i) random sampling of the Gaussian model of each state, and (ii) random sampling from the measured data. In this work, we choose the latter, random sampling from the measured data, since the Gaussian distribution may not accurately describe the randomness of the device's states and device to device variations.

B. MVM USING RCAs
RRAM crossbar arrays can perform the MVM operation, which is equivalent to n 2 multiply and accumulate (MAC) operations, with O(1) time complexity compared to O(n 2 ). The weight matrix is programmed/stored in the RRAM array cells as conductance values, and the input is applied as a voltage at the rows of the array. By grounding the columns of the array, the output current per column is proportional to the inner product between the input voltage vector and the conductance vector of the column, which can be written as where I j is the current of the j th column (i.e. post-synaptic current), G ij the synaptic weight in conductance, V i the i th input voltage (i.e. pre-synaptic voltage) and V n+1 = b which is the bias. The conductance of RRAM can only realize a positive value; however, both negative and positive weight realizations are needed in any neural network. In order to create negative weights, two weight realization techniques have been introduced: (i) using two RRAM cells per weight [32] as shown in Fig. 4, which is referred to as balanced realization, and (ii) using one RRAM as weight, in addition to one shared reference RRAM with the conductance of G r = (G max + G min )/2 ≈ G max /2, which is referred to as unbalanced realization [26], [33] where G max and G min are the minimum and maximum achievable conductances, respectively. In this work, we consider the first realization, which has double the dynamic range (conductance range ≈ (−G max , G max )), making it less susceptible to noise and variability at the expense of doubling the area and power. We would like to highlight that the realization, shown in Fig. 4, including analog partial sum and sum & compare circuits, is one way for realizing the MVM partitioning in the analog domain. Other works, such as ISAAC [6], use ADCs/ DACs to preforms the partial sum and sum & compare in the digital domain. We would like to emphasize that our framework is agnostic of the peripheral operations' realization since it focuses on IR drop problem. Besides, the IR drop has less impact on the overlap between the states, as will be discussed in Section III. The IR drop also causes each device conductance to be scaled differently, which could cause more dependency on the stored data. The differential output current can be written as where G ij is the differential conductance and can be written as matrix-vector multiplication as follows where V is the input voltage vector (e.g., input image) and the bias value. The current vector, I, is sensed and shaped by the activation function, which is mathematically described as where f (·) is the activation function. In practice, a DNN layer can be too large to be realized in hardware using a single crossbar array. Thus, these large layers are partitioned into smaller crossbar arrays connected to perform MVM as a single layer [6], [26]. In this work, we partition each layer into differential 128 × 128 arrays, which is the same size as recently fabricated arrays [34]. Larger array sizes lead to worse IR drops, causing higher degradation in performance, as discussed in [26].

C. WEIGHT MAPPING
Each weight is translated into a pair of conductance values, which can be mathematically formulated as where W max is the maximum value of the weight. If it is required to realize W max , G + and G − are set to G max and G min , respectively, and G = G max − G min . The difference between the two conductance values is constant and proportional to the required weight value, and each conductance is constrained to be between G min and G max .
Using the aforementioned 4-bit device, it is possible to realize up to 5-bit weight when two devices per weight are used. Table 1 and Fig. 5 show the weight mapping from quantized weight to device's conductance. Mapping-I is designed to occupy the entire dynamic range of the devices, which results in less overlapping for smaller precision. On the other hand, Mapping-II is designed to use the closest possible states to the zero state (i.e., high conductance state) with equal spacing. This results in less power consumption with higher overlap, as shown later in the results section. There is a linear relation between the device conductance and the weight if chosen properly except for the 5-bit case because of the high gap between the low conductance state and the first high conductance state. Thus, in this work, we consider up to 4-bit quantized neural network.

III. INEVITABLE WIRE RESISTANCE PROBLEM IN CROSSBAR STRUCTURES
The wire resistance is inevitable in nanostructure crossbar arrays. It is expected that the wire resistance would reach around 90 for 5 nm feature size [35]. The wire resistance creates undesired IR voltage drops that accumulate across the columns in the array leading to unwanted paths between the input and output nodes. Thus, the columns are not grounded anymore, resulting in a highly distorted MVM result. These voltage drops are a function of the stored data and the wire resistance.
Due to the lack of the experimental data for 128 × 128 crossbar arrays, we incorporated the device model discussed in the previous section into a SPICE-like simulator [35]. This simulator accurately simulates the interconnect and devices parasitics three orders of magnitude faster than SPICE with no loss in accuracy. Then, this simulator is incorporated in our framework for generating the training masks and for validating the performance of the retrained networks with the proposed method, which will be discussed in Section IV. It is also worth mentioning that any device model can be incorporated in our framework and the device impact on the neural network performance can be evaluated accordingly.
Using the analysis in [35], it can be shown that the output current with wire resistance effect can be modeled as (see Appendix A for the proof) where g(·) is the IR drop nonideality function, G is the programmed conductance matrix . u is the input vector and G e the effective weight matrix which has the same size as the RCA, which is n×m. In the balanced realization case, the output current is I = (G + e − G − e )u where the IR drop behaviour is the same for both RCAs. Thus, the difference can be seen as constant. On the other hand, in case of the unbalanced realization, the output current is I = g([G + , G − ])u) which is a more complex relation. Thus, the IR drop has severe impact on the output current in the unbalanced realization case. Fig. 6a shows the normalized measured effective weights for differential crossbar array with 1 wire resistance for two crossbar arrays, 256 × 256 array and 128 × 128 array, populated with random data. The measured weights decrease exponentially across the diagonal. Increasing the crossbar array size increases the IR voltage drop across the array, creating more sneak paths. Fig. 6b also shows the histogram of random data with 1-bit ternary quantization (i.e −1, 0, 1). The histogram of the 256×256 array has a smaller mean value and a larger standard variation compared to the 128 × 128 array. In addition, Fig. 6c shows the effect of the wire resistance on the histogram of the measured data for 1 , 5 , and 10 wire resistance for the 128 × 128 array. The higher the wire resistance, the more severe the IR drop. Figures 6 (d) and (e) show histograms of the quantized states for 3-bit and 4-bit weights for the 128 × 128 case. Ideally, each state should be a narrow pulse and non-overlapped, but because of the IR drop problem, it becomes wider and overlapped.
Partitioning each layer into the small arrays is necessary for three main reasons; • IR drop problem: partitioning helps to reduce the effect of the problem, compromising the main benefit of RRAMs, which is the density.
• Driver nonideality: each crossbar is driven by a driver circuit or buffer. The loading of the driver circuit is the input resistance of the driven row creating a voltage divider with the output resistance of the driver circuit. For example, the worst-case occurs when all the devices within the same row have a low resistance state (LRS), and the output resistance of the driver is R o . Thus, the voltage delivered to the crossbar input is So, it is necessary to reduce the number of devices per row, N , to mitigate the driver's effect. Otherwise, it has to be taken into consideration during the training. We do not consider it in this work for two reasons 1) with a well-designed peripheral circuit, its effect can be eliminated or mitigated. And 2) the IR drop problem is the main cause of the performance degradation [26]. However, we study its effect on the performance to find the required output resistance value of the driver circuit for the designers in section V.
• Fabrication problem: It is less complex to fabricate small crossbar arrays with high reliability. In addition, it is recommended to use two separate crossbar arrays for positive and negative conductances to have symmetric IR drop behavior. The corresponding conductances are scaled by the same value, unlike using a single crossbar array for both positive and negative conductances where each conductance will be scaled with different values.
To have accurate inference results, it is needed to run the inference with SPICE simulation with all circuits included. A SPICE simulator is adopted where the weight matrices are partitioned into small crossbar arrays as discussed in [26], [36] and simulated using a transient simulation for different input samples. Fig. 7 shows the simulation time of matrix-vector multiplication using SPICE for a 256 × 256 array partitioned into smaller arrays and for different input samples. The SPICE simulation time, without the peripheral circuits such as neuronal sensing circuit and drivers, increases exponentially with increasing the crossbar array size and linearly with increasing the input samples. On the other hand, the same Fig.shows the numerical SPICE-equivalent simulator adopted from [35]. The numeric simulator runs 140× faster than SPICE for one input sample and 1000× faster than SPICE for ten successive input samples. It is worth highlighting that the numerical results of the MVM are the same as the SPICE results.
Although the numerical model runs orders of magnitude faster than SPICE simulations, it is better not to be included in the DNN framework. It would take considerable training time even for small networks such as the MNIST dataset. Thus, it is better to have solutions that can be used to describe fabricated hardware. In this work, we use the numerical or SPICE-like simulator, without loss of generality, as a reference due to the lack of the hardware.

IV. PROPOSED IR-QNN TRAINING AND INFERENCE
A. QNN TRAINING FRAMEWORK Due to differential weight realization, it is possible to realize 2 n + 1 states where n ∈ {1, 2, 3, 4} is the device's precision in bits. We use binarized activation function {−1, 1} for simple and fast communication between the crossbar layers and eliminate area-and power-expensive blocks such as ADCs and DACs.
IR-QNN framework is an extended version of BinaryNet [37] to support more weight states and a bipolar activation function. During training, real-valued weights are quantized through stochastic rounding to the equally distributed point sets within [−1, 1]. Activations are binarized into {−1, 1} and physically implemented with ±0.1 V . Multilayer perceptron (MLP) and convolutional neural network(CNN) models are used for MNIST and CIFAR10 dataset, respectively (see Table 4 and 5). Convolution filters are 4D tensors but reshaped to 2D matrices with the number of output channels as the leading dimension so that convolution can be performed by matrix multiplication.
The modifications to the BinaryNet framework, to include higher quantized states and the IR drop estimation technique in the forward and backward computation, are shown in Algorithm 1. In the forward computation section, the modifications are as follows: (1) partitioning the quantized weights of each layer, W b k , into small weight matrices P b k , (2) application of the IR drop estimation method ( IRestimate function) to each partitioned array to create W e to be used in (3) the forward inference and (4) backward computation instead of the quantized weights W q . Similar modifications can be added to other QNN frameworks to capture the effect of the sneak path problem.
Accumulating the parameters gradients: ← λη end for

B. SYSTEM LEVEL ESTIMATION OF IR DROP PROBLEM
In this section, we introduce different methods to estimate the IR drop problem without the need for SPICE or numerical simulations.

1) TRAINING WITH MULTIPLICATIVE NOISE
It is clear in Fig. 6 that each state is spread with a certain statistical distribution due to the IR drops. One way to overcome this problem and to enable quick estimation of realistic weights, is to create a statistical model for each state and include it into the DNN framework which can be done as follows: where W e is the wire-resistance-effect-compensated weight matrix, W q i is the quantized weight matrix having i th state (such that i W q i = W q equals the quantized weight matrix), Q the number of states, element-wise multiplication and n i is multiplicative noise of i th state.
Each state is statistically modeled to different statistical distributions. The log-normal distribution is found to the best distribution (e.g., the highest likelihood) to describe the variability per state. The positive and negative states have similar histograms. Table 2 shows the curve-fitted model parameters for different wire resistance simulating different IR drop scenarios. Although the method would have the same effect on the summed current per column, this method is not very effective since it treats all the locations in the array equally, which does not describe the real behavior of the IR drop problem, as shown in Fig. 6a.

2) TRAINING WITH MASKS
Another solution is to generate average mask that can account for the cell location in the array. This mask is element-wise multiplied by the quantized weight matrix similar to [26], as follows: where M is the average mask matrix. This mask method helps to predict more realistic behavior of a crossbar array and can be easily calculated with fabricated crossbar array. The mask matrix is generated from either SPICE simulations or equivalent numerical methods [35] and is normalized to the ideal desired current. Masks are generated by averaging the results of many (e.g., 1000) SPICE simulations using random input weight matrices [26], [38].
Due to the averaging, the generated mask is static and fixed. However, in practice, W e has some stochastic behavior around this average mask. Thus, an additive white Gaussian noise can be added to the mask to exhibits more practical behavior, which we refer to as a stochastic mask.
The third solution is to generate a mask for each state and element-wise multiplied by its corresponding state matrix. This solution can be mathematically formulated as follows: where M i is the corresponding mask matrix of i th state. Fig. 8 shows mask examples for training a 2-bit neural network having 5 states per weight. Fig. 8a shows the M ±0.5 and M ±1 masks for 1 wire resistance as an example. Furthermore,, Fig. 8b shows the effect of the wire resistance on the M ±1 mask. It is clear the high degradation with higher wire resistance values. It is worth mentioning that applying masks during training has a negligible effect on the training time. In order to test the efficacy of estimating the IR drop effect, we calculated the mean square error (MSE) between the estimated effective weight matrix and the measured effective weight matrix for the single mask and multiple mask sets for 128 × 128 RCA with 10 . Clearly, the multiple mask set shows more than 100× better MSE compared to the single mask set, which is also is verified with DNN experiments later.

C. BATCH NORMALIZATION DURING INFERENCE
One of the important practices in training deep neural networks is adding a batch normalization layers to speed up the training, improve the performance and overcome vanishing gradient problem in DNNs, without the need for small learning rates which slow down the convergence [39]. In other words, the batch normalization is data whitening (i.e., removing the mean and variance of the data), which can be mathematically defined as follows for j th neuron where µ j and σ j are the mean and the standard deviation values of the input vector y j , and γ j and β j are trainable parameters. VOLUME 8, 2020 During the inference, these batch normalization layers can be removed by merging them with the preceding layers. As aforementioned, we use the sign activation function, thus the output activation o can be computed as Using sign(AB) = sign(sign(A)B) property yields Since σ always has positive value, batch normalized layer merged with the proceeding layer can be written as . The matrix form of (13) can be written as whereW = WD sign(γ ) and matrix whose diagonal is v. It is worth highlighting that the weight parameters are kept quantized even after merging batch normalization.

D. EXPERIMENTAL SETUP
To evaluate the effectiveness of our proposed technique we use the MNIST and CIFAR10 datasets [40], [41]. For each dataset, we use the same network architecture as given in BinaryNet [37], with these two changes: (1) instead of binary weights, we use up to 4-bit (or 17-state) quantized weights, (2) the model sizes are reduced. For MNIST the number of hidden neurons is reduced to one-fourth, and for CIFAR10, the number of channels is reduced to half, roughly reducing the number of model parameters to about 1/16 and 1/4, respectively. Furthermore, we also present results for modified AlexNet [42] for CIFAR10 dataset, with 3x channel sizes to justify our work on larger networks. Details of the networks can be found in Table 4, 5, and 6. The experiments are performed in three main steps. The first step is the baseline training, which uses floating-point weights/activations to obtain the best test accuracy, where test accuracy is the ratio of the correctly recognized samples for unseen data. We use the default training parameters for 100 epochs of MNIST training and 500 epochs of CIFAR10. Learning rates halve every 20 or 50 epochs for MNIST and CIFAR10, respectively, initially from 2 −6 . The baseline accuracies are 98.4% for MNIST dataset on multilayer perceptron network, 88.5% for CIFAR10 dataset on modified   VGG9 network and 81% for CIFAR10 dataset on modified AlexNet network, which is similar to the best accuracies for the datasets reported in the literature. The accuracy reported is test accuracy, that is, the inference accuracy for unseen data.
The second step is fine-tuning, which is running additional training iterations using the weight from the first step as the initial weight. While the weight in the first step gives a very high accuracy on GPU, it is unlikely to give good results if used for RRAM crossbar arrays due to distorted MVM computation caused by IR drop. The fine-tuning retrains the networks to mitigate the discrepancy by using the mask methods for different wire resistances and quantization levels. During fine-tuning, we use learning rates starting from 2 −9 , training additional 50/200 epochs for MNIST/CIFAR10 models. The other parameters remain the same as the baseline training. At the beginning of fine-tuning, the accuracy plummets due to the introduction of the mask but eventually recovers through fine-tuning. Note that the accuracy at the end of fine-tuning is not indicative of the real performance of RRAM crossbar arrays since it is not trained with SPICE simulations, for which we need a separate validation step.   The third step is validation. The output of the second step is quantized weights after merging the batch normalization layers to be programmed to RRAM crossbar arrays. The trained weights are mapped to the RCAs, as discussed in section II. The unused portions of the RCAs are populated with random data so that the distribution of IR drops do not change [26]. The required number of RCAs to implement each network is shown in Table 4, 5, and 6. For CNN, there are two ways to implement conventional layers either by flattening the convolutional layers into one large matrixvector multiplication, which results in the lowest latency at the expense of the area and power or by iterating over the same kernel multiple times, which results in high latency but saves area. In our framework, we use generalized matrix multiplication (GEMM) with im2col mapping for accelerating the convolutions layers where the filters and input patches are laid out into a 2-D matrices, which are multiplied to compute the same dot product of the convolutional operation [43].
Our validation setup takes the quantized weights, and runs SPICE-based RRAM crossbar simulation, to get the effective output currents with the device model that has been presented in Section II-A. The effective output currents are fed back to our QNN inference framework to obtain networklevel inference results. Note that neither training nor mask is used during validation. The test accuracy obtained from validation is what we can expect to see if the quantized weights are perfectly programmed to RRAM crossbar arrays, barring stochastic and other unmodeled noise/faults during RRAM read. Some of those nonidealities are considered in our additional experiments (see Section V).

V. RESULTS AND DISCUSSION
In this work, we consider three test scenarios, with 1 , 5 , and 10 wire resistances to consider different technology nodes. For instance, for a 50 nm feature size, it is expected to have around 5 wire resistance [24], [35]. The results shown in Table 7 illustrates that without considering the IR drop problem in training, the accuracy drops to around 10 ∼ 12% from the baseline test accuracies regardless of the number of bits.
After the retraining using the proposed techniques, the networks were able to reach close to the baseline accuracies. Tables 8-11 show the validation accuracy for MNIST and CIFAR10 datasets for different wire resistance scenarios and different mappings. The higher the wire resistance, the higher the drop in performance. Clearly, mapping-II shows a better performance for 1-bit and 2-bit cases since the used resistance values are much higher than the one used for mapping-I, and the severity of the IR drop problem is determined by the ratio between the device's LRS and the wire resistance value. The smaller the ratio, the more severe the IR drop problem. In general, training with multiple mask set achieves the best performance among the proposed solutions. On the other hand, increasing the number of bits does not improve the performance monotonically. The accuracy drops with increasing the number of bits to more than 2 bits, which is attributed to the overlap between states, as illustrated in Fig. 6. The drop in CNN test accuracy is much higher than MLP test accuracy due to its high sensitivity to weight variations. Clearly, from these results, the proposed training method provides the best performance with 2-bit devices. VOLUME 8, 2020   In Fig. 9, we compare our method against noise-injection Adaption (NIA) method [27] which was proposed for QNNs. The comparison is performed on VGGNet on CIFAR10 dataset with mapping-I. NIA method performance highly drops with increasing the wire resistance and with increasing the weight precision as well. Clearly, our method outperforms NIA method in all scenarios thanks to capturing the spatial behavior of the IR drop problem.

A. STUCK-AT FAULT EFFECT
Stuck-At Fault (SAF) defects cause another inevitable problem that affects the accuracy results of the MVM, which is the main operation in DNNs. SAF defects vary based on the fabrication technology and RRAM switching materials. In some recent works, the percentage of SAF fabricated crossbar arrays is about 10% for 1024 × 1024 for an in-house test array [44], [45] and is about 0.2% for 128×64 array (with only 15 devices stuck off) [46]. With the knowledge of the exact locations of the SAF devices, the network can be trained to isolate the SAF devices or at least mitigate their effect. However, this is not practical for DNNs where the trained weights are not designed for specific hardware. It should run without any knowledge of the location of SAF defects. Thus, in this work, we explore the effect of different SAF percentages on the recognition accuracy without retraining the network, assuming that the SAF devices are randomly distributed in each crossbar array. Fig. 10 shows the recognition accuracy with changing the SAF percentage for a 10 wire resistance scenario. The performance has no significant drop up to 50% and 20% stuckat open (OFF) for MNIST and CIFAR10 datasets, respectively. On the other hand, performance is more sensitive to stuck-at close (ON) case, where the stuck at close can cause weight-sign flip, a scenario that significantly affects performance and does not happen for the stuck-at open. For instance, for 1-bit, the possible weight values are {−1, 0, 1} which are mapped to G + = {HRS, HRS, LRS} and G − = {LRS, HRS, HRS}, respectively. The stuck at OFF would occur to one of the 2 LRSs, which would result in mapping all original weights to zero weight value in the worst case. On the other hand, the stuck at ON can occur to one of the four HRSs, which causes that mapping {1, −1} to ''0'' and ''0'' to be mapped to either ''−1'' or ''1''. Thus, the performance is more sensitive to the stuck-at ON. However, The effect of stuck at On can be mitigated by implementing ''0'' by two LRSs, which would hurt the stuck-at OFF performance. Thus, In our implementation, we choose two HRSs to implement ''0'' since stuck-at OFF is the most common in these devices [44]- [46]. The results in Fig. 10 shows that 1-bit and 2-bit networks are more robust against Stuck at ON compared to 3-bits and 4-bits cases. In addition, the stuck-at OFF results show some accuracy drop regardless of weight precision. The higher weight precision has a slight accuracy improvement after knee point. In conclusion, there is no need to retrain the network with the full knowledge of the SAF devices' locations, especially that the reported SAF percentages are less than 20% [44], [45].

B. EFFECT OF DEVICE VARIABILITY
In order to achieve low variation in each conductance's state, i.e., precise programming, write error-correcting techniques such as write-verify technique are usually used. However, such techniques require multiple write and read cycles, increasing the programming time of the entire crossbar array. In addition, a write-disturb problem occurs where writing some cells might disturb the written data in other cells [47]. With multiple writes to the same cell, a higher rate of the disturbed cell can occur. Thus, it is better to write once and take the variability into consideration during the training and validation. In order to consider these device variabilities, we have added Gaussian noise to each conductance's state where we vary the normalized standard deviation, σ n , normalized to the difference between the states ( g = 6µS) from zero to 100%. Fig. 11 shows the effect of changing the normalized standard deviation on the MLP and modified VGG networks using MNIST and CIFAR10 datasets, respectively. In this experiment, we have used the trained model with the multiple mask set and sampled from the measured data shown in Fig. 3a without performing any retraining to include new noise. Clearly, the performance slightly drops with increasing σ n . In case of the MLP network, the performance drops around 1.5% in the worst case for mapping-I. Approximately, the performance is the same regardless of the number of bits for mapping-I. On the other hand, mapping-II is more sensitive to the variations. The performance drops around 4% for 1-bit and 2-bit cases with 1 wire resistance after σ n = 0.5. This drop is caused by the narrow spacing between the states in mapping-II. Training with higher wire resistance values reduces the drop in the performance to 1.5%. In the case of a modified VGG network, the accuracy is slightly dropped by 2.1% at worst case compared to the zero noise case. Thus, network sensitivity to the programming variability is small since the IR drop problem dominates the programming noise model.

C. EFFECT OF LIMITED RETENTION
The main factor that would affect the performance of the DNNs over time (i.e., aging) is the RRAMs' retention. Recent works show different retention values based on device structure and materials. These works also show that the device drifts over the time towards a very low conductance state, G f , which is less than the formed low conductance state, G min [48]- [50]. In addition, the drift speed is a function of temperature. The higher the temperature, the higher the drift that occurs [48]. Due to the lack of aging modeled, we adopt the following retention model to be incorporated in the validation simulations to study the performance degradation with time. We emphasize that aging was taken into consideration in the training process to simulate practical scenarios. The conductance change versus time can be modeled as follows where G i is the initial conductance state, v d is the drift coefficient, and t n normalized retention time which is normalized to the retention value of the device. We chose this normalized model to simulate different RRAMs' behaviors with different drift coefficients. Fig. 12 shows the effect of aging on the validation accuracy of the MNIST dataset for two scenarios v d = 10 and 0.1 with 25% variability in the normalized retention time to simulate different device conditions. Clearly, the network was able to achieve the baseline accuracy for more than the 50% of the retention time. Then, we start seeing performance degradation regardless of the number of bits. The accuracy degrades faster for a smaller drift coefficient.

D. POWER AND AREA RESULTS
The power dissipation during the inference consists of the power dissipation of RCAs and the power dissipation in the peripheral circuits. However, due to the resistive nature of RCAs, the power is mainly consumed inside the RCAs. Fig. 13 shows the power dissipation of RCAs for processing one input image at 0.1 V . Clearly, mapping-II consumes VOLUME 8, 2020  around 20% and 35% of the power consumed in mapping-I for 1-bit and 2-bit cases, respectively. It is worth to highlight that the power consumption of using the model trained with the average mask is the same as the one trained with multiple mask set. To estimate the performance metrics, we used the hardware setup shown in Fig. 4, discussed in [36] in which direct communication between the RCAs without network on chip for routing purposes (i.e., fully dedicated hardware) which would give the highest performance. Adding a network on chip offers a highly flexible design; however, it costs power, area, and latency. Thus, with hardware setup, the total power consumption per image is estimated to be around 0.9 W, 132 W, and 225 W for the MLP, Modified VGG, and Modified AlexNet networks, respectively, where RCAs consumes 65.8%, 69.65%, and 91.5% of the total power while the rest is consumed in the sensing circuits and memory cells needed to store intermediate stages while pipelining. On the other hand, the total area is estimated to be around 0.0185 mm 2

VI. DRIVER AND NEURONAL CIRCUITS REQUIREMENTS
As previously discussed in section IV, we trained the networks to have binary activation function, {−1, 1}, for efficient communication and buffering between the fully connected layers. Three nonidealities need to be considered while designing the periphery circuits: • Driver output resistance: Each crossbar array is driven by a driver or buffer circuit. The output resistance of the driver circuit creates a voltage divider with the parallel RRAMs within the same driven row.
• Neuronal input resistance: After the current is summed within the crossbar array, a current sensing circuit is needed to sense the summed current from the positive array and compare it with the summed current from the negative array, and give a positive or negative output voltage. The input resistance of the sensing circuit creates extra loading to the crossbar array.
• Neuronal circuit variability: Due to PVT variations of the circuit, the comparison between current sensed from positive and negative RCAs is biased to one of them with a random value.
These nonidealities disturb the MVM computation, which affects the DNN performance. With well-designed circuits, there is no need to consider them during training. Ideally, the driver circuit should have zero output resistance, and the neuronal circuit should have zero input resistance and zero current offset.
In this section, we study the effect of these nonidealities on the MNIST network performance to find the maximum values that the network can tolerate without affecting the performance. It is worth to highlight that there is no retraining with peripheral circuits nonidealities. Including them in training  will relax the design requirements. Although, we see it is more beneficial to consider the worst case. Fig. 14 shows performance degradation due to the output resistance of the driver circuit with changing the number of bits and wire resistance. The results show that the trained network with multiple mask set can tolerate higher driver resistance, especially with a higher number of bits. In addition, mapping-II exhibits better behavior compared to mapping-I. Similarly, the effect of the load resistance is studied and shown in Fig. 15 for different wire resistance, number of bits, and different training methods. The performance degradation due to increasing the load resistance is the same for any number of bits. Training with a higher wire resistance value exhibits better performance against load and driver resistances. The effect of changing the standard deviation of the offset current on recognition accuracy is depicted in Fig. 16. Per these results, the network can tolerate up to 0.1 mA, 0.5 mA, and 1 mA standard deviation of the offset current for 1 , 5 , and 10 wire resistance, respectively, regardless of the number of bits per device.

VII. CONCLUSION AND FUTURE WORK
The paper proposed a QNN framework with a software-level technique to incorporate the IR drop problem in training the deep quantized neural networks. The efficiency of the proposed method is proven with three neural networks and is compared with prior work showing a significant improvement. We also studied the impact of other nonidealities, such as device variability, stuck-at faults, and aging. Our results show that the 2-bit device exhibits the best performance and training with multiple mask set. Our experimental results recommend using a driver with output resistance to be less than 200 , and input resistance of the sensing circuit should be less than 100 . In addition, the input-referred current offset deviation should be less than 0.1mA.
The main limitation of the proposed method is that the masks are generated assuming random data stored in RCAs which generate a static mask which might show lower performance than expected for nonrandom data patterns. Adding random noise to the average mask improves the performance in some cases, mainly for MLP networks and for Convnets with 1-bit weight precision case only, as shown in our experiments. In addition, the proposed method shows less performance with increasing the wire resistance for instance less the accuracy dropped around 5%p at 10 wire resistance. In this work, a software evaluation of the framework's performance is performed, taking into consideration the hardware limits. A full circuit validation with a full implementation of the peripheral circuits, such as partial sum, sample & compare, and max-pooling circuits, is also needed to validate the performance. Besides, other datasets should be tested in our framework. These two points are left for future work. Figure 17 show the crossbar array with wire and capacitive parasitics. According to [35], the nodal voltages can be obtained by solving the following 1 st order system of differential equations:

APPENDIX A STEADY-STATE MODEL OF MVM USING CROSSBAR ARRAY
where V and u are the nodal voltage and the excitation vectors, respectively and M, N and G are the coefficient matrices containing the RRAM's conductances, wire and capacitive parasitics values. The construction of these matrices can be found in detail in [35]. Since our concern is the steady-state behaviour, the capacitive parasitics can be ignored. Thus, the steady-state nodal voltage vector can be written as The output current is needed to have accurate MVM as discussed. Thus, the steady-state output current equation can be defined as I ss = V ss (18) VOLUME 8, 2020 where I o is the output current vector and is the selection matrix which is given as where R sBL is the parasitic load resistance of the crossbar array and m and n are the array dimensions. Consequently, the output current can be written as This equation can be written as I ss = G e u where G e = M −1 G Similar analysis can be adapted for nonlinear switching devices.  Dr. Eltawil has been on the technical program committees and steering committees for numerous workshops, symposia, and conferences in the areas of low power computing and wireless communication system design. He received several awards, as well as distinguished grants, including the NSF CAREER grant supporting his research in low power systems.