A fast algorithm to simulate nonlinear resistive networks

Authors: benjamin  scellier
Venue: NeurIPS_cc_2024_Workshop_MLNCP
Published: 2024-09-02 12:49:18

Abstract: Analog electrical networks are explored as energy-efficient platforms for machine learning. In particular, resistor networks have recently gained attention for their ability to learn using local rules such as equilibrium propagation. However, simulating these networks has been challenging due to reliance on slow circuit simulators like SPICE. Assuming ideal circuit elements, we introduce a fast simulation approach for nonlinear resistive networks, framing the problem of computing its steady state as a quadratic programming (QP) problem with linear inequality constraints. Our algorithm significantly outperforms prior approaches, enabling the training of networks 327 times larger at speeds 160 times faster.

Main paper: 
 introduction: 
 As energy costs for machine learning continue to rise, neuromorphic computing platforms are gaining attention as alternatives to traditional neural networks and GPUs (Momeni et al., 2024) . These platforms use analog physics and compute-in-memory architectures to achieve substantial energy savings. Among them, nonlinear resistive networks have sparked particular interest (Kendall et al., 2020; Dillavou et al., 2022 Dillavou et al., , 2024;; Wycoff et al., 2022; Anisetti et al., 2024; Stern et al., 2022 Stern et al., , 2024;; Kiraz et al., 2022; Watfa et al., 2023; Oh et al., 2023) . These networks use variable resistors (e.g., memristors) as trainable weights, combined with diodes for nonlinearity and voltage/current sources as inputs. Nonlinear resistive networks can be trained by gradient descent via equilibrium propagation (Kendall et al., 2020) or variants (Dillavou et al., 2022; Anisetti et al., 2024) and they are universal function approximators (Scellier & Mishra, 2023) . They leverage the principles of electrical circuits (e.g., Kirchhoff's and Ohm's laws) to perform inference and extract the weight gradients. Furthermore, the learning rules governing conductance changes are local. These features make nonlinear resistive networks good candidates as power-efficient learning-capable hardware, with recent experiments on memristive networks suggesting a potential 10,000x gain in energy efficiency compared to neural networks trained on GPUs (Yi et al., 2023) . Small-scale networks based on these principles have been successfully built and trained on datasets like Iris (Dillavou et al., 2022) , validating the approach.Larger-scale simulations of nonlinear resistive networks on tasks such as MNIST (Kendall et al., 2020) and Fashion-MNIST (Watfa et al., 2023) further highlight their potential. However, these simulations, often conducted using the general-purpose SPICE circuit simulator (Keiter, 2014; Vogt et al., 2020) , are hindered by SPICE's slowness. For example, Kendall et al. (2020) reported that simulating a one-hidden-layer network with 100 hidden nodes on MNIST took one week for just ten epochs of training. Due to the lack of efficient methods for simulating nonlinear resistive networks, some researchers have resorted to linear networks (Stern et al., 2022 (Stern et al., , 2024;; Wycoff et al., 2022) , thus missing an essential feature of machine learning: nonlinearity.To address this, we introduce a novel methodology for simulating nonlinear resistive networks. Our algorithm, suitable for networks with arbitrary topologies, is an exact coordinate descent algorithm for convex quadratic programming (QP) problems with linear constraints (Wright, 2015) . Applied to the 'deep resistive network' architecture of Kendall et al. (2020) , our algorithm becomes an 'exact block coordinate descent' algorithm, whose run-time on GPUs is orders of magnitude faster than SPICE. The key contributions of this manuscript are:• We show that, in a nonlinear resistive network, under an assumption of ideality of the circuit elements, the steady state configuration of node electrical potentials is the solution of a convex minimization problem: specifically, a convex quadratic programming (QP) problem with linear inequality constraints (Theorem 1).• Using the QP formulation, we derive an algorithm to compute the steady state of an ideal nonlinear resistive network (Theorem 2). Our algorithm, which is an instance of an 'exact coordinate descent' algorithm, is applicable to networks with arbitrary topologies.• For a specific class of nonlinear resistive networks called 'deep resistive networks' (DRNs) (Kendall et al., 2020) , we derive a specialized algorithm to compute the steady state (Section 3). It exploits the bipartite structure of the DRN to perform exact block coordinate descent, where half of the coordinates are updated in parallel at each step of the minimization process. Each step of our algorithm involves solely tensor multiplications, divisions, and clipping, making it ideal to execute on parallel computers such as GPUs.• We perform simulations of DRNs, trained with equilibrium propagation (EP) on the MNIST dataset (Section 4). Compared to the SPICE-based simulations of Kendall et al. (2020) , our DRNs have up to 327 times more parameters (variable resistors) and the training time per epoch is 160 shorter, resulting in a 50000x larger network-size-to-epoch-duration ratio. We also train DRNs of two and three hidden layers. 
nonlinear resistive networks: 
 We study electrical circuits composed of voltage sources, current sources, linear resistors, and diodes, collectively referred to as nonlinear resistive networks. These networks are modeled as graphs, with branches containing individual elements. We denote the set of branchesB = B VS ∪ B CS ∪ B R ∪ B D ,to represent voltage sources, current sources, resistors, and diodes, respectively. For each branch (j, k), the voltage across a voltage source is v VS jk , the current through a current source is i CS jk , and the conductance of a resistor is g jk .We assume that these circuit elements are ideal, with behavior determined by the current-voltage (i-v) characteristics in Figure 1 . The network's steady state is defined by branch voltages and currents satisfying branch equations and Kirchhoff's laws. Under the assumption of ideality of the circuit elements, the steady state of a nonlinear resistive network is described as follows. Theorem 1 (Convex QP formulation). Consider a nonlinear resistive network with N nodes, and denote v = (v 1 , v 2 , . . . , v N ) the vector of node electrical potentials. Under the assumption of ideality, the steady state configuration of node electrical potentials, denoted v ⋆ , satisfiesEQUATIONwhere E : R N → R and S are defined as:EQUATIONS := {(v 1 , v 2 , . . . , v N ) ∈ R N , v j ≤ v k ∀(j, k) ∈ B D , v j = v k + v VS jk ∀(j, k) ∈ B VS }. (3)Some comments are in order. First, S is the set of feasible configurations of node electrical potentials, and E is the power dissipation function (sometimes also referred to as 'energy function' although this terminology is not physically accurate). Importantly, not all feasible configurations satisfy the laws of electrical circuit theory; Theorem 1 states that among all feasible configurations, the one that is physically realized (the steady state) is the configuration that minimizes E.Second, the feasible set S is defined by linear constraints. The inequality v j ≥ v k ensures nonnegative voltage across diodes (v j > v k if the diode is in the off-state, and v j = v k if the diode is in the on-state), and the equality v j = v k + v VS jk ensures fixed voltage across voltage sources. Certain conditions on network topology must be met to ensure S is non-empty. For instance, a loop of voltage sources with non-zero net voltage would violate Kirchhoff's Voltage Law (KVL), rendering S empty.Third, the function E represents half the power dissipated in resistors, plus the power dissipated in the current sources. Theorem 1 extends the 'principle of minimum dissipated power', stating that in a network of ideal voltage sources and linear resistors, the physically realized configuration minimizes power dissipation in the resistors.Finally, E(v) is a convex function of the node electrical potentials v, being a quadratic form with a positive definite Hessian. The feasible set S is also convex, defined by linear constraints. Thus, the steady state is the solution to a convex quadratic programming (QP) problem.Voltage Sourcei - + v0 v i v0 Linear Resistor i g v v i i = gv Diode i v v i Current Source i0 v v i i0Figure 1 : Ideal circuit elements and their current-voltage (i-v) characteristics. A linear resistor follows Ohm's law: i = gv, where g is the conductance (g = 1/r, with r being the resistance). An ideal diode operates in two states: "off" (i = 0 for v ≤ 0) and "on" (v = 0 for i > 0). An ideal voltage source maintains a constant voltage v 0 , independent of the current i, while an ideal current source maintains a constant current i 0 , independent of the voltage v.v 1 v 2 v 3 v 4 v 5 v 6 v 7 v 8 v 9 v VS 13 v VS 12 g 25 g 58 i CS 57 - + - +Figure 2 : A nonlinear resistive network. The voltage sources form a tree (in blue), so if we set e.g. v 1 = 0, we can immediately infer v 2 = -v VS 12 and v 3 = v VS 13 . The steady state is computed using exact coordinate descent (Theorem 2) on internal node electrical potentials (in black). For example, at node k = 5, p 5 is calculated from connected resistors and current sources, p 5 = (g 25 v 2 + g 58 v 8 + i CS 57 )/(g 25 + g 58 ), and v 5 is determined by clipping p 5 within diode-imposed bounds, v 5 = max(v 3 , min(p 5 , v 4 )). This process repeats for all nodes until convergence.Using the QP formulation (Theorem 1), we derive a numerical method to compute the steady state v ⋆ of an ideal nonlinear resistive network. Starting from an initial configuration v (0) ∈ S, our algorithm iteratively minimizes the function E(v) by generating a sequence of configurations v (1) , v (2) , . . . within the feasible set S, ensuring that E(v (t+1) ) ≤ E(v (t) ) at each step.Each step of the algorithm proceeds as follows. Let v = (v 1 , . . . , v N ) ∈ S. First, we pick an internal node k of the network. The function E of Eq. ( 2) is a quadratic function of v k with other node electrical potentials fixed. Specifically, E takes the formE = a k v 2 k + b k v k + c k , where a k := 1 2 j∈BR g jk and b k := -j∈BR g kj v j -j∈BCS i CS jk .The function E(v k ) is bounded below and its minimum in R is obtained atEQUATION)The range of feasible values for v k , however, is constrained by the diodes. It is of the formEQUATIONTheorem 2 (Exact coordinate descent). The configurationEQUATIONIn particular v ′ ∈ S and E(v ′ ) ≤ E(v).Using Theorem 2, we minimize E via an 'exact coordinate descent' strategy. Starting from a feasible configuration v ∈ S, we iteratively select an internal node k, compute the value v ′ k that minimizes E, and update the configuration. This process is repeated until the energy E converges to a minimum. 
deep resistive networks: 
 The deep resistive network (DRN) model, introduced in Kendall et al. (2020) , is a layered nonlinear resistive network inspired by neural networks. DRNs hold promise for hardware implementation with memristive crossbar arrays. In this section, we highlight the advantages of DRNs for simulations on multi-processor systems like GPUs.In a DRN, a set of nodes connected to ground by voltage sources forms the input layer, while another set of nodes forms the output layer, whose electrical potentials represent model outputs.Power dissipation function. DRNs are inspired by the architecture of multi-layer perceptrons (Figure 3 ). Each node in the network, also called a 'unit', has an electrical potential, analogous to a unit's activation in a neural network. We denote v (ℓ) k the electrical potential of the k-th node of layer ℓ. Nodes in consecutive layers are interconnected by variable resistors, denoted g(ℓ)jk , serving as trainable weights. The power dissipation function of a DRN is:EQUATION)Feasible set. Nonlinearity in DRNs is introduced by placing diodes between unit nodes and ground, creating excitatory (non-negative potential) and inhibitory (non-positive potential) units. In hidden layers (1 ≤ ℓ ≤ L -1), even-indexed units are excitatory, and odd-indexed units are inhibitory. The output units (ℓ = L) are linear, without diodes. The feasible set is:S = {v ∈ R L ℓ=1 N ℓ | v (ℓ) k ≥ 0 if k is even, v (ℓ) k ≤ 0 if k is odd, 1 ≤ ℓ ≤ L -1, 1 ≤ k ≤ N ℓ }. (8)Input voltage sources. A key challenge in resistive networks is the non-negativity of weights: a conductance is non-negative. This is partially addressed by using both excitatory and inhibitory units in hidden layers. To further enhance representational capacity of the DRN, the number of input units is doubled to N 0 = 2 dim(x), where dim(x) is the input dimension. The input voltage sources are set such that vEQUATION2k-1 for each 1 ≤ k ≤ N 0 . Furthermore, to prevent voltage decay as network depth increases, input voltages are amplified by a factor A ≫ 1, setting v 0 2k-1 = +Ax k and v 0 2k = -Ax k .Output switches. For equilibrium propagation (EP) learning, each output node is connected to ground via a switch, a resistor, and a voltage source used to set target values. These output resistors share a common conductance value β > 0. During inference, switches are open; during training, they close to drive output units toward target values.In a DRN, the update rules (Theorem 2) for the 'hidden' units take the following form:EQUATIONAs for the output units, the update rule is:EQUATIONwhere β = 0 if the output switches are open.To speed up the computation of the DRN's steady state, we can exploit the fact that each v(ℓ)k update depends only on layers ℓ -1 and ℓ + 1. Thus, all units in layer ℓ can be updated simultaneously. We obtain an 'exact block coordinate descent' algorithm. We can push this idea further, by partitioning the network into even-indexed and odd-indexed layers, and updating these groups alternately.Equations ( 9) and ( 11) can be expressed in matrix-vector form, making each step of the algorithm involve ∼ L matrix-vector multiplications, ∼ L/2 divisions, and ∼ L/2 clipping operations. This makes the algorithm highly efficient for parallel computing platforms like GPUs. 
simulations: 
 We used our exact block coordinate descent algorithm to train deep resistive networks (DRNs) with equilibrium propagation (EP) (Scellier & Bengio, 2017; Kendall et al., 2020) on the MNIST classification task. We trained DRNs with one, two, and three hidden layers (1024 units each), labeled DRN-1H, DRN-2H, and DRN-3H. Additionally, we trained two single-hidden-layer DRNs: the DRN-XS model (100 units) of Kendall et al. (2020) , and the DRN-XL model (32,784 units), which is 327 times larger. Compared to the SPICE-based DRN-XS simulations of Kendall et al. (2020) , the DRN-XL model was trained for 10 times as many epochs (100 vs 10), 16 times faster (10 hours and 27 minutes vs one week). Thus, our network-size-to-epoch-duration ratio is 50000x larger. The DRN-XL model also achieves a significantly better test error rate (1.33% vs 3.43%).We also trained DRNs with backpropagation (BP) as a baseline. Although EP approaches the performance of BP, the small gap can likely be attributed to EP only approximating the gradient of the cost function, rather than extracting the exact gradient like BP does. Experimental details are available in the full version of the manuscript (Scellier, 2024) , and the code for reproducing the results is available at https://github.com/rain-neuromorphics/ energy-based-learning. 
discussion and conclusion: 
 Although this study focuses on layered architectures (DRNs), our algorithm applies to arbitrary network topologies, including unstructured (disordered) networks (Stern et al., 2022 (Stern et al., , 2024;; Wycoff et al., 2022) . Parallelization in such networks is also possible, though less straightforward to implement. Additionally, while we used equilibrium propagation (EP) for training, our exact coordinate descent algorithm is compatible with other learning methods, such as 'coupled learning' (Dillavou et al., 2022) , 'frequency propagation' (Anisetti et al., 2024) and 'agnostic EP' (Scellier et al., 2022) . More generally, by framing the steady state as a quadratic programming (QP) problem with linear constraints, our approach is compatible with various other QP-solving methods like primal-dual interior point methods (IPM) and sequential quadratic programming (SQP).Our methodology to simulated nonlinear resistive networks presents limitations too. To simplify realworld complexities, it assumes ideal circuit elements (resistors, diodes, voltage sources, and current sources). Real diodes, however, exhibit non-ideal behaviors such as forward voltage drop and reverse leakage current. Despite these simplifications, our algorithm can be valuable for prototyping and exploring various network topologies. Future work could extend our approach to simulate networks with more realistic i-v characteristics, such as piecewise linear curves. (However, it remains unclear whether smooth nonlinear characteristics such as exponentials could be efficiently and accurately simulated with this approach). Another caveat of the proposed DRN architecture is the need to amplify input voltages to prevent weight decay across the layers. This amplification can be very large, e.g. A = 4000 for the 3-hidden-layer DRN -see Scellier (2024) for the details of experiments. However, this can be mitigated by using 'bidirectional amplifiers' at each layer, instead of one large input gain (Kendall et al., 2020) .Looking ahead, our simulation methodology could significantly accelerate research in nonlinear resistive networks. It enables large-scale simulations of deep resistive networks (DRNs), providing a pathway to assess the scalability of these analog self-learning machines on more complex tasks. 
 
 Appendix: Second Workshop on Machine Learning with New Compute Paradigms at NeurIPS 2024 (MLNCP 2024). 



Reviews: review: title: title: A useful tool for simulating analog ML networks
main_review: main_review: This paper describes a means for simulating nonlinear resistive networks for performing machine learning in a way that vastly outpaces SPICE simulations. The results are clearly shown and the algorithm is well-described. It is easy to see this being useful for determining the ideal working properties of a prospective analog machine learning circuit. 


I have only a few minor comments and questions:

An indication of the required ranges of conductance and of voltage (max/min) would be helpful to assess the feasibility of some the architectures mentioned. 

Piecewise linear functions are mentioned in the discussion – are smooth nonlinear functions such as exponentials (a more accurate model of a diode) only feasible with this algorithm if approximated piecewise?

Would an optimizer like Adam be useful for speeding this optimization process? Do the authors select initial conditions they know to be close to the minimum (e.g. the solution to the previous training step)?
title: title: Paper suggests a new way of efficiently solving the responses of certain non-linear resistor circuits, providing for efficient computational implementation of learning algorithms in such systems.
main_review: main_review: Electronic circuits with adaptive resistor elements emerged recently in theory and experiment as a possible approach to implementing fast and energy efficient neuromorphic computing. However, computationally efficient ways of simulating such networks are lacking.
This paper suggests a quadratic programming approach (with linear inequality constraints) to efficiently compute the responses of these networks, containing certain ideal circuit elements including linear resistors and diodes. The addition of diodes as inequality constraints break the linearity of the network and allows it to compute, in principle, any nonlinear function. The authors explain how quadratic programming is useful in simulating such networks, implement coding tools for a certain topology reminiscent of deep neural networks, and demonstrate large speedups over previous computational approaches involving detailed circuit simulations (e.g. Spice).

This is overall an excellent paper that significantly extends the computational ability to simulate resistor-based learning circuits, and I can clearly see why this would be very useful for upcoming research on this topic. The ideas are clearly articulated and justified, showing why QP is natural for solving this problem. The simulated results also clearly indicate the computational advantage of this approach over previous ones. I also liked the authors idea of optimizing over the weights of odd and even indexed layers simultaneously. I am happy to recommend this paper for the MLNCP workshop.

I have some comments for the authors' consideration below:
1) Eq. 2 and following text: This is not an energy function but a power dissipation function, which is minimized due to Onsager’s law. While in theory this system behaves like an energy-based model, it actually isn’t really based on energy minimization. This is somewhat important as it means the network requires constant energy input to operate, as opposed to an actual thermodynamic “energy-based” system.
2) For the general coordinate descent algorithm, would you randomly pick an internal node potential at each iteration? Could one design a “smart” way of picking these nodes to accelerate convergence, i.e. by estimating which nodes require large changes to locally minimize power?
3) Is input amplification necessary in DRNs? If resistance values are small, I’d expect only mild voltage decay across layers.
4) Table and preceding text: Was it expected that backpropagation would give different results compared to equilibrium propagation? Both compute the same gradient and should scale similarly in terms of runtime. Also, please indicate units for wall clock time.