Stochastic gradient descent (SGD) is a widely-used, first-order optimization method for training neural networks. During forward pass, a neural network infers output and estimates loss by comparing the output with expected values. Based on the loss, gradients of weights in a neural network are calculated during backward pass and later updated during the weight update phase by error back-propagation algorithm. The amount of weight update is proportional to the calculated gradient with the scaling factor called learning rate, η:

where w_ij indicates a weight from pre-synaptic neuron i to post-synaptic neuron j, L is the loss of a neural network, and ∇_{ij, k}L indicates the derivative of L with respect to w_ij. x_i is the input activation of pre-synaptic neuron i, and δ_j is the error of post-synaptic neuron j propagated from output neurons.

Analog computation with a resistive cross-point array resembles fixed point or quantized number system unlike the digital computation done with floating point number system. For the forward pass, the input activation, x_i is converted to a corresponding discrete voltage, xiq, through analog-to-digital converter (ADC). Then, x_i is encoded into a pulse, PWM(xiq), by pulse width modulation (PWM). As the pulse trains are applied to the rows, currents flowing through the devices are summed at each column line, and one can obtain the weighted sum, z_j:

where operation “*" includes the integration of currents through the time domain and g_ij is conductance value of i^th column and j^th row device. Then, analog-valued z_j is again quantized into zjq by ADC, and zjq propagates to the next layer as inputs. For the backward pass, the error of pre-synaptic neuron i, e_i, is calculated in a similar fashion to the errors of post-synaptic neurons, δjq :

In the update phase, xiq is converted to a pulse train, SPG(xiq), by stochastic pulse generator (SPG) for stochastic update of a resistive memory device. In compliance with the rules in Equation (1), the update of g_ij is determined by x_i and δ_j.

where SPG function is defined by pulse bit length (BL) and parameters for pulse probability, satisfying η = η_xη_δ. η_x and η_δ can be modulated, such that two SPGs will have similar pulse generating probability (Gokmen and Vlasov, 2016; Gokmen et al., 2017). Δg_ij(g_ij) is the amount of conductance change of g_ij as a function of g_ij, which is dependent on the update characteristics of the specific memory device. Here, we assume a linear model for Δg_ij(g_ij) which is described in section 2.2.

When mapping a neural network into resistive cross-point arrays for analog operations, we note that each weight can be expressed with a single or multiple memory cells depending on the architecture or weight mapping scheme (Xiao et al., 2020). As an example, a pair of memory cells, g_ij and g_{ij, ref}, can be read differentially to represent a single weight: w_ij = K(g_ij−g_{ij, ref}), where K is a factor that scales the weight into the conductance (Gokmen and Haensch, 2020).

Practical resistive memory devices feature various types of non-ideal characteristics such as asymmetric conductance response (or update asymmetry), short retention time, device-to-device variation and limited number of states. Among them, update asymmetry is known to hamper the convergence (Huang et al., 2020) of SGD-based neural network training and causes a significant performance drop. As shown in Equation (4), it is the Δg_ij(g_ij) term which distorts the actual update from the ideal update behavior defined in Equation (1). That is, g_ij will represent a weight value distorted form ideal software calculation after several cycles of weight update unless Δg_ij(g_ij) is unity. For most variants of resistive memory devices, as illustrated in Figure 1A, Δg_ij(g_ij) scales up or down the amounts of updates, which results in inconsistent updates of the conductance of devices (or weights). This type of behavior not only misguides the direction of updates, but also reduces the effective number of states in resistive memory devices.

Figure 1A illustrates the conductance(g) and the conductance change behavior (Δg) which can be described by the following equations, Equation (5):

g_{min, ij} and g_{max, ij} is the minimum and maximum conductance state of resistive device and Δg0,ij± is the amount of Δg_ij when g_ij of the device is at g_{0, ij} = (g_{max, ij}+g_{min, ij})/2. s_ij, namely slope, is defined as Δg0,ij±/gmax(min),ij for nonlinear device, which represents the reciprocal of approximate number of states of the device. For ideal device, i.e., perfectly symmetric and linear device, slope is defined zero. d_ij,k is number of pulse needed to achieve desired amount of w_ij,k (see Supplementary Materials section 1.1 for the details).

To embody the update asymmetry of a device, we introduce a parameter called asymmetry factor, AF(>0), which is associated with existing parameters: Δg0,ij=AFij·Δg^0,ij, where Δg^0,ij is a unit amount of Δg_{0, ij}. Thus, Equation (6) can be rewritten into Equation (7):

Equation (7) shows that symmetry of device is solely defined with AF·Δg^0,ij, if g_{min, ij}, g_{max, ij} and Δg^0,ij± is fixed. Therefore, we can compare the asymmetry of device within the same range by modulating AF. In the case of ideal device, slope equals zero and Δg_{0, ij} is constant regardless of AF_ij. However, to make the notation consistent, we denote this case as AF_ij = 0. For instance, there are three different cases in Figure 1A, namely, AF_ij = 0;AF_ij = 2.0;andAF_ij = 6.0. All three cases have the same g_{min, ij}, g_{max, ij} and Δg^0,ij. We note that many devices show non-zero AF behavior (Brivio et al., 2018; van De Burgt et al., 2018; Islam et al., 2019; Kim et al., 2019b).

Tiki-Taka algorithm (Gokmen and Haensch, 2020) is proposed to resolve the performance degradation issue in neural network training caused by the update asymmetry in the crosspoint elements. We showed that the asymmetric update prevents devices from being accurately updated by the amount of calculated gradients. Figure 1B illustrates the key difference between SGD and Tiki-Taka algorithm in the weight update phase. In the case of SGD, an array, W, stores the weight vectors of a neural network. In the update phase, a weight, w_ij, is updated with the gradient ∇_ijL(= x_iδ_j). On the other hand, Tiki-Taka algorithm requires one additional array, namely A, and it stores ΔW by accumulating gradient vectors, ∇L. The weight vectors stored in the array A are denoted as W^A and the array, C, stores the weight vectors denoted as W^C. When compared to SGD, Tiki-Taka algorithm accumulates the gradient, ∇_ijL, in wijA. Then, the accumulated gradient is transferred to update the weight, wijC at every ns steps. Tiki-Taka algorithm supports sparse update by adopting cell-selection vector, u_π, where u_π(i) is the i^th element of a sparse vector, u, such as a one-hot vector and π is a permutation, π:{1, ⋯ , n} → {1, ⋯ , n}. Finally, Tiki-Taka algorithm allows the array A to participate in the forward pass by introducing a parameter, γ∈[0, 1]: wij=γwijA+wijC. Therefore, the effective weight w_ij depends on γ: if γ = 1, wijA+wijC replace the effective weight vector. Otherwise, if γ = 0, then the effective weight is wij=wijC. In summary, Tiki-Taka algorithm can be formulated as follows:

where ∇ijγL indicates the derivative of the L with respect to wij=γwijA+wijC. η and λ indicates SGD learning rate and transfer learning rate, respectively.

Here, we formulate the impact of asymmetry by adopting even-odd function decomposition, with which any function can be represented (Gokmen and Haensch, 2020), and integrate the model into the SGD update rule in Equation (1) (see Supplementary Materials section 1.2 for the details):

where the function Sym represents symmetry of the devices and the function Asym corresponds to asymmetry. The weight, w_ij, is scaled with a factor K from the conductance, g_ij: w_ij = K·(g_ij−g_{ref, ij}) and Δw^0,ij=K·Δg^0,ij. w_{max, ij} is a correspond weight to the maximum amount of conductance, g_{max, ij}. To analyze it in the aspect of optimization, Equation (10) is reformulated as follows:

Definition 3.1 (Mismatch factor (MF)). Mismatch factor (MF) is an additional factor appeared in the original SGD weight update rule, originated from the update asymmetry of the weight storage device:

MF is multiplied to the gradient of a weight, ∇_ijL, and affects the weight optimization process. As observed in the Definition 3.1, MF_ij scales with AF, and is dependent on the sign of gradient and the current weight state, w_ij, of a device. The dependence on the sign of gradient can create fluctuation of MF_ij if sgn(∇_ijL) changes the sign, and AF_ij magnifies the fluctuation of MF_ij. Therefore, there can be abrupt change in MF_ij affecting the neural network training.

Figure 2 displays the classification error in MNIST handwritten digit recognition problem as a function of training epoch for various scenarios of training algorithm and AF values. We perform the experimental simulation using the open-source toolkit aihwkit (Rasch et al., 2021) and further details about the experimental environment can be found in Supplementary Materials section 3.3. First, performing SGD with asymmetric devices results in severe accuracy drop both for the train and test datasets when compared with the software baseline. Since the update asymmetry leads to non-unity MF and consequent unfavorable behavior in weight update process, a sharp drop in performance is found during the training. On the other hand, the cases with Tiki-Taka algorithm and finite asymmetry show robust performance reaching to the floating point baseline (Gokmen and Haensch, 2020). As shown in Figure 2, when AF^A = AF^C = 1.0, the train and test error gradually decrease and can reach to the baseline. However, the remaining question is how much asymmetry Tiki-Taka algorithm can tolerate. Other cases shown in Figure 2 with various pairs of AF^A and AF^C answer the question. The train error does not gradually decrease, which indicates that the weight vectors of neural networks are stuck at bad local optima. Interestingly, the case with AF^A = AF^C = 0 shows higher train and test error than those of SGD with asymmetric non-linear devices. This result indicates that there exist optimal pairs of AF, and a reliable method to determine the optimal pair of AF^A and AF^C is necessary. To find the method, it is necessary to quantify and understand how MF affects the optimizer, Tiki-Taka algorithm, during the neural network training.

Recent variants of SGD introduce the concept of history-dependent update and utilize gradient information at each step, including momentum and adaptive gradient, to overcome the issues of vanilla SGD (Kandel et al., 2020). In such methods, the weight update is determined not only by the current gradient calculation, but also by the history of previous gradients and their importance at the current step. For instance, if the current weight update includes only the first-order gradient, then the history dependent update is composed of first-order gradient vectors and their respective importance for each step.

Definition 4.1 (Gradient Scheduler). History-dependent optimization algorithm is generalized by the governing rules including n^th order gradients:

t = t(k) (k = 1, ⋯ , T) is a gradient scheduler or scheduler, which contains information on the importance of the k^th gradients, sgn(∇ijL),∇ijL, ⋯,(∇ijL)n.

Therefore, the scheduler, t(k), defines how the gradients of previous steps for k < T affect the update of current step at k = T, since t(k) contains information on weighting of the gradients dependent on the step, k.

From the Definition 4.1, we can find the t(k) for Tiki-Taka algorithm as follows (see Supplementary Materials section 2.2 for the detailed derivation):

where 1_T≡ns0 is an indicator function which returns 1 every ns steps, t(k) is controlled by three hyper-parameters (Gokmen and Haensch, 2020): γ, ns, and u_π. These parameters cause the weight matrix W to be updated sparsely and asynchronously.

Figure 3 compares the gradient schedulers in three algorithms: vanilla SGD, momentum-based optimizer, and Tiki-Taka algorithm. The parameters of Tiki-Taka used are ns = 1, u = 1, and γ = 0, and the ideal, symmetric and linear, switching characteristics are assumed for the devices in both arrays (i.e., AFijA=AFijC=0). With ideal devices, Tiki-Taka algorithm memorizes all the previous gradients with equal importance. Therefore, it resembles the momentum-based optimizer if β → 1. However, it should not be confused with the momentum-based optimizer, because the scheduler of Tiki-Taka algorithm with asymmetric devices depends not only on the step, but also on the history of gradients and slopes of devices unlike the momentum-based optimizer.

By including the update asymmetry, we can reformulate the scheduler of Tiki-Taka algorithm, Equation (14), as follows (see Supplementary Materials section 2.2) for the detailed derivation):

Equation (15) indicates that both MF^A and MF^C play an critical role in weight optimization. MF^A kicks in every step of k, and MF^C provides a factor until the last step of T in the scheduler of Tiki-Taka algorithm. Therefore, the optimizer of Tiki-Taka algorithm is uniquely determined by two MF's. Considering that wij,TA is determined by the history of gradients, ∇Lij,1γ,⋯,∇Lij,Tγ, the scheduler of Tiki-Taka algorithm is parameterized by the step and the history of the first-order gradients. MF^C affects the optimization process for different T, T+1, ⋯  steps in common by scaling all the ideal t(k). Therefore, MF^C is one of key components in building the optimizer, and MF^A determines the optimizer of Tiki-Taka algorithm.

As we discussed in the previous sections, the update asymmetry in array A and C plays an important role in neural network training with Tiki-Taka algorithm as an optimizer. To experimentally explore the impact of AF^A and AF^C on training neural networks, we perform a series of experiments by varying AF^A and AF^C values and training a multi-layer perceptron (MLP) with MNIST dataset. For this experiment, we set (w_min, w_max) = (–1,1) and Δw^0=0.001 for all devices in the array as mentioned in section 2.2, which illustrates the device reaching w_max = 1 from initial weight value of w_min = −1 within 600 updates if (xiq, δjq) = (1, –1) is given during update phase (see Supplementary Materials section 3.2 for the details of experimental environments).

Figure 7 shows the heat map of train (Figure 7A) and test accuracies of MLP as a function of AF^A and AF^C values and with two different set of SGD and transfer learning rates (Figures 7B,C). The results show the consistent observations as analyzed in section 5. First, when AF^C is fixed, the train and test accuracies peak at AF^A = 1.78 and decrease as it becomes larger or smaller. Smaller AF^A causes the optimizer to accumulate the history of gradient with equal importance, which hampers the convergence of the weight vectors due to oscillation. However, this effect is reduced with larger AF^C thanks to the adaptive gradient impact. In contrary, larger AF^A values disallow accumulation of the full history of gradients, which results in reduction of the regularization effect. This tendency is more significant for larger AF^C values, where the variance and expectation of updates for the array C is large enough not to store the entire history of gradients. Considering that larger AF^A causes memorization of short-term gradients, the regularization effect becomes weakened.

The optimal combination of AF values, which provide the best accuracy, can be changed when the combination of learning rates vary. For example, in Figure 7B, the test accuracy peaks at the combinations, (AF^A, AF^C) ∈ {(1.0, 1.0), (1.78, 1.0), (1.0, 1.78)} when the transfer learning rate and SGD learning rate is (0.01, 0.02). However, as shown in Figure 7C, the combination of AF is optimal when (AF^A, AF^C)∈{(3.16, 1.0), (5.62, 1.0), (10.0, 1.0)} if the neural network is trained with different pair of learning rates, (0.02, 0.01). Considering that AF of each device is typically pre-determined at the moment of device fabrication, it is crucial to find the robust AF^A and AF^C values over the learning rate space which do not degrade the accuracy. From this perspective, finding the robust region of AF without searching the best pair of learning rates can be guaranteed by introducing the thresholds of certain measurements such as test accuracy. The following robustness score provides the hard-bound threshold for finding the optimal combination of AF over the space of learning rates.

Definition 6.1 (Robustness score). Suppose that m is a neural network model, and Meas(·) is a measurement such as accuracy or loss of the model m after training for a given dataset D. H is a hyper-parmeter space for software training, H={(η,λ), where η = SGD learning rate andλ = transfer learning rate}. A Robustness score, RS(m), for the given threshold, th, is defined as follows:

Figure 8 illustrates the robustness score, RS(m), over a learning rate space, H={(η,λ)}, where η∈{0.01, 0.02, 0.04} and λ∈{0.01, 0.02, 0.04}. As discussed in section 3.2, AF of each array affects the test accuracy of the network using Tiki-Taka algorithm. Test accuracies do not reach sufficient value when AF of each array is too small (AF^{A, C} = 0) or large (AF^{A, C} = 10). Therefore, optimal AF value pairs which maximize the robustness score exist near AF^{A, C} = 1 as shown in Figure 8. The score map has two implications. First of all, if the minimum requirement of test accuracy (threshold) is determined, the region whose robustness scores are over a certain criteria could be used to define the minimum specifications of update asymmetry. Therefore, one can obtain the range of AF value for each array which gurantees the minimum test accuracy. Second, by scanning over various pairs of learning rates, one can identify the region which provides the AF values which are most independent on the learing rates. To find such region, the score is post-processed with interpolation and gaussian filtering. The robustness score is approximated by piece-wise linear interpolation in the domain of AF^A and AF^C. The interpolated data is smoothed by gaussian filtering (g(AF^A, AF^C))) to obtain the robust region in the domain: g(AF^A, AF^C) = (2π^σ2)−1/2·exp(−(((A^FA)2+(A^FC)2)/(2σ²)), where we define the degree of “robustness" of the score map by modulating σ of a gaussian filter. Finding the robust region with a gaussian filter is consistent with the definition of variations in AF of devices. In Figure 8C, the robust region is around AF^A = 1.2 and AF^C = 1.0 with the score greater than 8.95. As we experiment with a total of 9 scenarios of learning rates, the region with the score greater than 8.95 displays the required test accuracy within the given range of learning rate pairs. When compared to the region with RS(m) = 4.5, the region with RS(m)≥8.95 shows almost doubled efficiency during on-device training. In other words, it requires half of the resources to find the optimal pair of learning rates. This calibration method to find the best combination of AF is applicable to other types of hyper-parameters if the hyper-parameters space H can be expanded to include them.

In practical applications, the target performance of the network, such as the test accuracy of the MNIST classification task using MLP, is one of the most significant design choices as some tasks require the maximum accuracy while others do not. Considering that the robustness score, RS(m), varies as target accuracy changes, studying the relationship between RS(m) and the target threshold value can provide insights on the choice of AF for each array. Figure 9 illustrates the robustness score maps with two different threshold values, 97.0 and 97.5%. When the threshold is set to 97.0%, the score is robust against the change in AF^C value and the efficiency of searching over the hyper-parameter space, H with AF^A is not degraded in the proper region. However, in the case of the threshold being 97.5%, the best choice is to fix AF^C around 1.0. Then the asymmetry requirement for the array A is lifted. Although the robustness scores of the case with the threshold being 97.5% are relatively lower than those of 97.0%, selecting AF^A and AF^C of the region has comparative advantages over the other regions.

Consistent with the analysis in sections 5 and 6.1, Figure 9A describes that the neural networks are able to find the optimum with a wide range of AF^C values over H if AF^A is in the proper region. This observation is due to the adaptive gradient effect that AF^C brings, which helps the convergence around the nearest local optimum. In this case, the neural network converges to the local optimum instead of further exploring the loss surface. This results in the high robustness score region along the wide range of AF^C, though it does not guarantee higher accuracy. In Figure 9B, on the other hand, the region satisfying over 97.5% accuracy is spread along the axis of AF^A, which indicates that the specification for AF^C value is strict. If AF^C value is within the specification, then the robustness score depends on the choice of AF^A. With smaller AF^A values, the optimizer allows the weights of neural networks to explore over loss surface, which enables to find better optimum.

Recent research on SNN and its implementations using resistive memory devices have shown that the devices are designed explicitly to be embedded with the learning algorithm for SNN such as Spike-timing-dependent plasticity (STDP) and equilibrium propagation (Scellier and Bengio, 2017). Accordingly, the update asymmetry of devices is more likely to affect the process of training neural networks (Kwon et al., 2020). Inspired by this motivation, they demonstrated that asymmetric non-linear devices are more powerful than symmetric linear devices (Brivio et al., 2018, 2021; Kim et al., 2021). Our analytical tools and calibration method can be applied to illustrate the relationship between the update asymmetry of devices and its impact on training in detail. Therefore, it is worthwhile to pursue further studies on the relationship in more extensive range of cases of learning algorithms.