Memristor is an umbrella term for describing a broad class of memory technologies that follow a state-dependent Ohm’s law (Chua, 2014). Physical realization of memristors comes in several forms, including resistive random access memory (ReRAM), spin transfer torque RAM, phase change memory, ferroelectric RAM, and others (Chen, 2016). In essence, these devices store a non-volatile conductance (or resistance), which can be modified by providing a large write voltage and can be read using a smaller read voltage. The conductance is bounded between 2 extreme values, G _min and G _max. Memristors are particularly attractive for neuromorphic computing because they exhibit behavioral similarity to biological synapses, combining storage, adaptation, and physical connectivity in one device. Moreover, combining multiple memristors into high-density crossbars enables the efficient computation of vector-matrix multiplication (VMM), where the (voltage) input vector to the crossbar columns are multiplied by the matrix of memristor conductances to produce the (current) output vector.

Huge numbers of VMM operations are performed in neural networks during training and inference. When implementing neural network weights as memristor conductances in hardware, there will be no need for sparse design off chip weight storage and data movement as in most of digital based designs (Jouppi et al., 2017; Davies et al., 2018). This yields high energy efficiency, which is an important factor for AI on edge devices (Lee and Wong, 2016). Computation on information based VMM can be represented by currents, voltages (voltage mode), or a combination of the two, each approach has its own set of strengths and weaknesses (Merkel and Kudithipudi, 2017). Voltage-mode based VMM circuits, are the most common approach, in which the inputs are represented as voltages, and the outputs are represented as currents. The current-mode VMM (Merkel, 2019) has some advantages such as low supply voltage, current-mode design techniques, etc., but current distribution can be challenging (Marinella et al., 2018; Sinangil et al., 2020). Charge-based VMM is another approach aiming to perform dot product operation, using voltage inputs and to charge binary-weighted capacitors and performing summations through charge redistribution among capacitors via switched capacitor circuit principles (Lee and Wong, 2016). The main advantage of this approach that there is no static power in the circuit and there is no limitation in technology node down scaling. However, multiple clock cycles are needed to perform the multiplication operation. Time-based VMM approach is another way to implement addition in the analog domain by using a chain of buffers. The delay of each buffer in each stage can be modified according to the weight input summation for each stage (Everson et al., 2018). In time-domain computing, the values are encoded as discrete arrival times of signal edges (Freye et al., 2022). Another implementation of time-based VMM is discussed in (Bavandpour et al., 2019; Sahay et al., 2020), in which memristor-dependent currents are summed and integrated on a capacitor and then the charge is converted back to the time-domain representation.

Quantization methods for deep learning are becoming popular for accelerating training, reducing model size, and mapping neural networks to specialized hardware. The simplest quantization methods use rounding to reduce activation and weight precision after training. This usually results in large drops in accuracy between the full-precision and quantized models. Other methods quantize weights, activations, and sometimes gradients during training, resulting in better performance (Hubara et al., 2017). In this work, we only quantize weights and activations. The core idea is to use quantized values during forward propagation and full-precision gradient estimates during backward propagation. For activations, we use a simple threshold model on the forward pass: x = 1 2 sign s + 1 2 (1) where sign (⋅) is 1 if the argument is non-negative and −1 otherwise. Since the sign has a gradient that is zero everywhere it will stall the backpropagation algorithm and nothing will be learned. To fix this, we approximate the gradient as ∂ x ∂ s ≈ 1 1 + exp − k s 1 − 1 1 + exp − k s , (2) where k was empirically chosen as 2. In other words, on the backward pass, the gradient is calculated as if the activation had been a logistic sigmoid function. Of course, we note that the threshold activation function is indeed a logistic sigmoid with a k value of +∞.

For weights, we use the following quantization technique: w q = 2 × r o u n d 2 Q − 1 c l i p w , − 1,1 + 1 2 2 Q − 1 − 1 (3) where Q is the desired number of bits for the weight, round (⋅) rounds to the nearest integer and clip (w, a, b) = max (a, min (b, w)), where a ≤ b. For backpropagation, we estimate the gradient as ∂J/∂w ≈ ∂J/∂w _q

The core of the design uses domino logic style neuron based on memristor RC delays. In essence, memristors are used as a configurable RC delay. The delay of the memristor RC circuit represents synaptic weight computations and a simple binary neuron activation function represented by an inverter, as shown in Figure 1A. The operation of the synaptic weight matrix is divided into two phases: pre-charge and evaluation. When the clock signal ϕ is low, the dynamic node v _d (input to the inverter) is pre-charged to V _dd through a PMOS transistor. When the clock is high, evaluation starts and the dynamic node discharges at a rate depending on the pull-down network (RC time constant). Once the node reaches the threshold value of the inverter, the neuron’s output will go high. During the evaluation phase, the voltage on the dynamic node of neuron i in layer l evolves as v d i l t = V d d × exp − ∫ 0 t G i l ξ C d d ξ (4) where G i l is the equivalent pull down conductance. A memristor only contributes to the pull-down conductance when its selector transistor is on. Assuming that memristor conductance values are constant during the evaluation phase and input voltages are binary values, i.e., v x j l − 1 ∈ 0 , V d d , then G i l is a piecewise constant function written as: G i l t = 1 V d d ∑ j = 1 N l − 1 v x j l − 1 t G i j l (5) where N _l−1 is the number of neurons in the previous layer, plus 1 to account for the bias input.

We use a simple parasitic model for the capacitance in (Eq. 4), where each 1T1R synapse contributes one unit of capacitance C to the dynamic node from the NMOS drain. Assuming the PMOS transistor has minimal sizing, and the inverter has 2:1 PMOS:NMOS sizes ration, the total capacitance is estimated as C d = 4 + N l − 1 C = 4 + N l − 1 A m i n k o x ϵ 0 t o x , (6) where A _min is the minimum transistor area, k _ox is the relative permittivity of SiO₂, ϵ ₀ is the permittivity of free space, and t _ox is the transistor gate oxide thickness. Bounds on the time that it takes to discharge the dynamic node to the inverter threshold will be important to set an appropriate clock frequency. From (Eqs. 4–6), the minimum discharge time will be t m i n l = − ln θ V d d 4 + N l − 1 C N l − 1 G m i n = β N l − 1 G m i n (7) where θ is the threhold voltage of the inverter.

Notice that this expression is approximately independent of N _l−1 as N _l−1 becomes large, meaning that these neurons are self-normalizing. That is, the maximum excitability of a neuron from the activity of one-presynptic neuron is constant regardless of fan-in. In contrast to the minimum discharge time, the maximum time will be infinite if we ignore leakage through selector transistors. However, we need to ensure that the dynamic node does not take longer to discharge than the evaluation period. Otherwise, there will be information lost. Therefore, we set a bound on the bias, which corresponds to a 1T1R cell that has a constant gate voltage of V _dd : G i 0 l ≥ β T / 2 , (8) where T is the clock period. This ensures that t m a x l = T / 2 .

Figure 2 shows how the leakage current could affect the neuron behavior with increasing numbers of synapses. Here, we assume all of the synapse transistors are off, and each of the conductances is set to the maximum value, which will maximize the leakage current. The plot shows the final voltage on the dynamic node during a 50 ns evaluation period. As the number of synaptic input increases, both the capacitance on the dynamic node and the total leakage current increase. This has the effect of maintaining a relatively large voltage on the dynamic node, even with a large synaptic fan-in. In fact, as the number of synapses becomes large, the increased capacitance dominates, causing the effect of leakage to decrease. In all cases tested, the leakage current was not enough to discharge the dynamic node to the inverter’s threshold.

To represent positive and negative weights, two of domino logic style neurons are used, inhibitory neuron to represent negative weight components and excitatory neuron to represent positive weight components as shown in Figure 1B. The time difference between nodes (inhibitory and excitatory) to reach the threshold of the inverter can be given as: Δ t i l = t e x i l − t i n i l = β V d d 1 v x l − 1 ⋅ G e x i l − 1 v x l − 1 ⋅ G i n i l (9) This time difference encodes the input to the neuron’s binary activation function. Figure 1C shows an example where the excitatory domino circuit discharges faster than the inhibitory, leading to a positive value of Δt.

In this paper, we are interested in binary neurons, so it will be sufficient to know if Δ t i l is negative or positive, which can be calculated using an arbiter circuit, the details of which are discussed later. First, though, it is important to point out how (Eq. 9) corresponds to pre-synaptic neuron outputs, weights, and the post-synaptic neuron inputs. The neuron input is Δ t i l , and we define the unitless version of it as s i l = Δ t i l / ( T / 2 ) ∈ [ − 1,1 ] . The vector of pre-synaptic neuron outputs are v x l − 1 , and the unitless version is x l − 1 = v x l − 1 / V d d ∈ [ 0,1 ] . Finally, the weights are related to the conductances G e x i l and G i n i l . There are infinite ways to map a weight w i j l to two conductances G e x i j l and G i n i j l . In this work, we use a power optimized scheme, where minimum conductance values are used for the excitatory component and inhibitory component when the weights are negative and positive, respectively. Then, a linear function maps the opposite component: G e x i j l = G m i n + G m a x − G m i n m a x 0 , w i j l (10) G i n i j l = G m i n − G m a x − G m i n m i n 0 , w i j l (11)

Now, (Eq. 9) can be rewritten completely in terms of unitless values as: s i l = β ′ 1 x l − 1 ⋅ w e x i l − 1 x l − 1 ⋅ w i n i l (12) where w e x i l = G e x i l / G m a x and w i n i l = G i n i l / G m a x , and β′ = 2β/(TG _max ). Contrast this with the usual dot product input of a neural network: s i l = x l − 1 ⋅ w i l (13)

In other words, the neuron input of our design has a non-linear relationship with pre-synaptic neuron outputs and the weights, whereas typical formulations of neural networks have a linear relationship. This non-linearity stems from the natural non-linearity of RC circuits, and to remove it would require an additional evaluation phase, such as the method proposed in Bavandpour et al. (2019). Instead, we opt to keep the hardware as simple as possible and we note that the non-linearity could be accounted for in two ways. On one hand, the behavior of (Eq. 12) could be modeled directly in Tensorflow. While this would be the simplest solution, we have found that this leads to several simulation challenges. For example, it is easy for (Eq. 12) to have undefined values when inputs or weights become small. In addition, we observed unstable learning and, in some cases, inability to converge when working directly with (Eq. 12).

Interestingly, though, because our neurons are binary, only the sign of s i l is important, and we can use (Eq. 13) for training. Figure 3 compares the values in (Eq. 9), (Eq. 12), and (Eq. 13). Here, we performed 1000 Monte Carlo simulations with uniformly-distributed weights between −1 and 1, and Bernoulli-distributed inputs with probability value 0.5. From the plot, one can observe that the normal dot product operation in (Eq. 13), which we refer to as “software” maps non-linearly to the excitatory-inhibitory time difference as well as the normalized version of the time difference, referred to as “hardware”. The ranges of the hardware values are inversely proportional to the ranges of the software values, which is expected due to the inverse relationships in (Eq. 9) and (Eq. 12). Critically, all of the datapoints are in the lower-left and upper-right quadrants, meaning that the sign of the software and hardware data are always the same. In other words, the non-linear behavior of the neuron’s input will not affect its output. The next section discusses how this sign is captured using an arbiter to yield a binary activation function.

We explored two possible designs for implementing the arbiter-based activation function that converts the time difference between the excitatory and inhibitory domino circuits into a binary value: v x i l = 0 , Δ t i l ≤ 0 V d d , Δ t i l > 0 (14) Initially, we designed the arbiter as shown in Figure 4A. Here, NOR gates are used to disable the inhibitory (excitatory) circuit from discharging once the excitatory (inhibitory) circuit crosses the inverter threshold. An advantage of this approach is that only one of the domino circuits will fully discharge. For example, if the excitatory domino circuit discharges more quickly, it will cause the inhibitory circuit to go back into the pre-charge phase before it reaches the inverter threshold. This will reduce dynamic power consumed during the pre-charge phase. However, we observed that this approach has poor stability (seen in the waveform in Figure 4A especially for small Δt, and often led to oscillations as well as incorrect outputs.

The second design is shown in Figure 4B, which uses cross-coupled NAND gates and includes a metastability filter to ensure the output doesn’t remain long in an invalid logic state and consequently, it is possible for that design to have that behaviour and enter metastability. In addition, NOR gate two large PMOS transistors in series, which means more capacitance and more delay compared to the NAND gate. The main advantage of this design is that the memristor domino circuits are not included in the feedback path, so the circuit’s stability is not dependent on the memristor states. This is the design that we used for the rest of this paper.

In order to implement large multi-layer neural networks information from one layer needs to be transferred to the next layer. Synchronizing the transfer of information between layers is critical, and here we explore three techniques with various tradeoffs. In this section, we use the XOR problem as a case study, where our multi-layer perceptron (MLP) neural network consists of 2 inputs, 2 hidden neurons, and 1 output. The output should be ‘0’ when both inputs are the same and ‘1’ when the inputs are different.

Simulation of large neural networks in electronic solvers like SPICE is prohibitively slow. We use a combination of SPICE (Synopsys HSPICE) and Tensorflow in order to have fast simulations while still capturing key hardware behavior. Our simulation strategy is shown in Figure 11. The key aspects of the circuit behavior discussed in the previous sections have been captured using HSPICE simulations with a predictive technology 130 nm bulk CMOS transistor model (https://ptm.asu.edu/) and memristor parameters based on the memristor proposed in (Prezioso et al., 2015), which has ON and OFF conductance values of G _min = 5 × 10⁻⁷ and G _max = 5 × 10⁻⁵ and programming voltage magnitudes near 1 V. In this work, only a subset of the device conductance range is used: G _min = 1 × 10⁻⁶ S and G _max = 1 × 10⁻⁵ S. Since the device has approximately linear behavior at low voltages that are below the programming threshold, we have modeled it as a resistor. Tensorflow is used to train the neural network to get the weights, which are converted into conductances based on the weight mapping in (Eq. 10) and (Eq. 11). Finally, another Tensorflow model with identical neural network topology but more accurate hardware behavior (e.g., stochastic behavior of neurons, conductance-based weights, etc.) is used to estimate the hardware performance on the dataset.

We tested the proposed design approach for an MLP with 1000 hidden neurons on the MNIST handwritten digit dataset (http://yann.lecun.com/exdb/mnist/). The results are shown in Figure 12A. Without considering noise, the software and hardware simulations produce almost identical results, with accuracies in the high 90s. These accuracies hold approximately constant across different levels of weight precision, from 3 to 10 bits. Note that 1- and 2-bit precision results were much lower and are not included in the plot. As noise effects are added to the simulation, the accuracy generally decreases. Low and moderate noise results are almost identical, while high noise gives a clear drop in accuracy. However, even with high noise magnitude, the degradation is less than 2%. We have also explored the effects of conductance variations on the proposed hardware design. These variations may arise from imprecise programming or drift of conductance values over time. Small conductance variations (e.g., 10%) have negligible effects on the test accuracy while larger variations (20% and more) can start to have considerable effects. Note that the effects of conductance variations are more pronounced at higher levels of weight precisions, which may motivate employing lower-precision devices or fewer of devices’ conductance levels.

The power consumption of the proposed design was modeled by assuming that most of the power is consumed when a neuron pre-charges. The justification for this is that, especially for neurons with high fan-in the switching capacitance of the neuron’s dynamic node will be much larger than the capacitance at other nodes in the circuit. Therefore, the power can be formulated as P ≈ 3 × 1 + η ∑ l = 2 L α C l V d d 2 f (16) where η is a fitting parameter that comes from the extra power associated with the inverter, arbiter, etc., α is the switching activity factor, L is the number of layers, and C ^l is the total switching capacitance of the layer. For the dynamic pipelining synchronization scheme, α = 1, since each layer will pre-charge every clock cycle. In addition, the value of C ^l is 3C times twice the number of synapses in a the layer (to account for both excitatory and inhibitory). The factor of 3 comes from each synapse’s source, drain, and memristor capacitance. We have empirically found η ≈ 0.19. In Figure 13, we show the power consumption for 100 randomly-sized 3-layer networks vs. the number of synapses and neurons in the network. For each network, both the inputs and weights were generated randomly. Furthermore, the network used a clock frequency of 10 MHz. From this data, we estimate the energy efficiency of our design to be approximately 1.26 fJ per classification per synapse. A comparison with similar works that designed MLPs for MNIST classification is shown in Table 1. Our work has slightly better energy efficiency than that reported in Yakopcic et al. (2015) while giving much better accuracy. 2T2R are needed to represent positive and negative weights. For our work, the number of transistors per neuron is 20. For Yakopcic et al. (2015) we estimate the number of transistors per neuron to be 4, and the number of transistors per neuron in (Jiang et al., 2018) is 30 (15 for excitatory and 15 for inhibitory).