A scalable method to find the shortest path in a graph with circuits of memristors
AA scalable method to find the shortest path in a graph with circuits of memristors
Alice Mizrahi,
1, 2, ∗ Thomas Marsh, Brian Hoskins, and M.D. Stiles National Institute of Standards and Technology, Gaithersburg, USA Maryland NanoCenter, University of Maryland, College Park, USA (Dated: September 14, 2018)Finding the shortest path in a graph has applications to a wide range of optimization problems.However, algorithmic methods scale with the size of the graph in terms of time and energy. Wepropose a method to solve the shortest path problem using circuits of nanodevices called memristorsand validate it on graphs of different sizes and topologies. It is both valid for an experimentallyderived memristor model and robust to device variability. The time and energy of the computationscale with the length of the shortest path rather than with the size of the graph, making this methodparticularly attractive for solving large graphs with small path lengths.
I. INTRODUCTION
Optimization problems, such as sequence of decisions,resource allocation, or navigation, are omnipresent andimportant. They can be solved by being mapped to find-ing the shortest path between two nodes in a graph. Al-gorithms to do this exist, but the time and energy theyconsume scale with the size of the graph and thus be-come prohibitive for large systems.
In this paper, wepropose a method to solve the shortest path problem thatscales consistently better than with the size of the graph.Our approach is inspired by biological systems that needto solve optimization problems to be energy efficient. Forexample, members of ant colonies work in parallel tofind the shortest path to their food without supervision. Related algorithms, called ant colony optimization, takeinspiration from this phenomenon to solve optimizationproblems in an approximate but efficient way.
Morebroadly, the field of swarm intelligence provides optimiza-tion algorithms inspired from animal populations.
However, these algorithms are limited when they run onconventional hardware, which computes in a mostly se-quential way.Pershin and Di Ventra proposed to solve the short-est path problem directly in hardware, using nanodevicescalled memristors, which have dynamics that provide re-inforcement mechanisms similar to the ones at play inthe ant colony optimization.
Memristors are definedby having conductances that change when subjected toelectrical current.
As a voltage is applied to a net-work of memristors forming a graph, more current willflow through the shortest branch because that branch hasthe lowest resistance. This increased flow will cause anincrease of conductance of the memristors on the short-est path, attracting even more current. When steadystate is reached, the memristors on the shortest pathhave a much larger conductance than the ones on thelonger paths, making it possible to electrically read outthe shortest path.In this work, we build on Pershin and Di Ventra’s ideaand propose a modified method. The advantage of thismethod is that it does not require prior knowledge aboutthe shortest path. We show that the time and energy consumed by this method scale with the length of theshortest path rather than with the size of the graph,which makes it potentially more efficient than algorith-mic methods.Sections II and III provide background on memristorsand describe how they can be used to find the short-est path in a graph. Sections IV and V present ourmethod and validate it through numerical simulations oflarge numbers of randomly generated graphs of varioussizes and topologies. We show that it is valid for realis-tic experimentally-derived memristor model and param-eters, and that it is robust to device variability. SectionVI addresses the key advantage of this method by show-ing that the time and energy it consumes scales withthe length of the shortest path. Finally, Section VII ad-dresses the hardware implementation of this method.
II. MEMRISTOR-BASED OPTIMIZATION
Memristors are a class of devices that exhibit hystereticbehavior: their electrical conductance can be modified ina non-volatile fashion. The conductance can take val-ues between two extreme states. The device is said tobe ”ON” or ”OFF” when in its highest or lowest con-ductance state, respectively. The conductance can berepeatedly increased and decreased by running currentthrough the device.
A graph can be represented by a circuit of memristors:the nodes of the graph are electrical junctions connectedby memristors implementing the edges of the graph. Anexample is shown in Figure 1(a): each black dot corre-sponds to a node and each colored rectangle to a mem-ristor. The goal is to find the shortest path between thestart and end nodes, marked by green stars in Figure1(a).A voltage is applied across these nodes, causing cur-rent to flow through the circuit. The system is initializedwith all memristors in their OFF state. The shortestpath, i.e., composed of the lowest number of memristors,is more conductive than the other paths. As a conse-quence of Kirchhoff laws, more current flows through thememristors on the shortest path than through the other a r X i v : . [ c s . ET ] S e p memristors. This causes their conductance to increase,drawing even more current to them and thus creating areinforcement mechanism. On the other hand, the mem-ristors outside the shortest path have little current flow-ing through them so their conductance does not increase.This behavior is shown in Figure 1(b) which presentsthe evolution of the conductance of the memristors be-longing to the shortest path (in red) and the others (inblack) versus time. After some time, the system reachesa steady state where the memristors on the shortest pathare ON while the others are OFF. This state is depictedin Figure 1(a), where the conductance of each memristoris represented by the color of the rectangle. The suc-cessfully found shortest path can be visually observed inred, contrasting with the rest of the graph in purple. Astheir state is non-volatile, the voltage across the circuitcan be turned off and the individual memristors can bemeasured in order to determine the shortest path (as de-scribed in Section IV).Note that the shortest path is found here without su-pervision: it emerges from the dynamical evolution of thesystem. Furthermore, the computation is parallel, witheach memristor evolving at the same time. This paral-lelism is key to how the time and energy consumed by thecomputation scale as the size of the problem increases,as detailed in Section VI. III. MEMRISTOR MODELS
In this work we consider both a simple generic memris-tor model as well as a more complex realistic memristormodel, derived from experiments. We start by studyinga simple generic memristor model, where I = V ( G ON x + G OFF (1 − x )) (1)where I is the current though a memristor and V thevoltage across it. G ON = 10 − S and G OFF = 10 − Sare the extreme conductance states of the memristors. x is an internal state variable, bounded between 0 and 1,describing the conductance dynamics of the device. Theconductance can be increased by running current throughthe device. dxdt = γ | I | − xτ (2)where γ = 10 A − s − , and τ = 0 . et al ., with the device pa-rameters they extracted from experiments. The model corresponds to a Pd / WO / W stack. Applying a voltageacross the device induces oxygen vacancy migration inthe oxide. The width and length of the created conduct-ing filament determine the conductance of the device,through the internal state variable x . I = (1 − x ) α (1 − exp ( − βV )) + xγ sinh ( δV ) (3) dxdt = λ (exp ( η V ) − exp ( − η V )) − xτ (4)where α = 5 × − S, β = 0 . − , γ = 4 × − S, δ = 2 V − , λ = 4 . − , η = 0 .
004 V − , η = 4 V − and τ = 10 s are device parameters.Note that in this model, the polarity of the voltagematters, as it can grow or shrink the filament. Since mostoptimization problems map to directed graphs, devicessensitive to the polarity of the voltage could be an asset.For broader applications, schemes where two memristorsof opposed polarities are connected in parallel for eachedge could be used, as proposed in Ref. 16. Furthermore,other types of memristors, such as phase-change memo-ries, are not sensitive to the polarity of the voltage acrossthem. The non-linear form of the current-voltage re-lationship complicates the solution of the system of equa-tions, leading us to use a commercial circuit simulator tosolve the system of equations.
IV. CONTROL VOLTAGE
We examine the importance of the control voltage, i.e.,the voltage applied across the circuit, and show that us-ing a constant voltage is not practical for applications.Panels (a), (c) and (d) of Figure 1 show the steady stateof one single graph after application of three differentvoltages, as described in Section II. The generic linearmemristor model described in Section III was used. InFigure 1(a), where a voltage of 4 mV was applied, thecomputation was successful. However, Figures 1(c) and1(d) show examples of incorrect computation. In Figure1(c) the voltage (1 mV) was too low to significantly in-crease the conductance of any memristor. On the otherhand, in Figure 1(d), the voltage (10 mV) was too high sothe conductance of some memristors on longer paths wasincreased comparatively to those on the shortest path.Only an optimal range of control voltage leads to a suc-cessful computation of the shortest path.To investigate the role of the control voltage on thecomputation, we define a metric of success for solvingthe shortest path problem. The input of the problemis the topology of the memristor circuit, i.e., the graph.In order to read the output of the computation (oncesteady state is reached or after the evolution of the sys-tem is stopped by the user, as we describe in Section V),we proceed as follows. The shortest path starts with thestart node – where the voltage source is applied – andis constructed node by node. The conductance of each G ON C ondu c t an c e Time (s) G OFF
Shortest path G ≈ G ON Other memristors G ≈ G OFF Δ G ≈ Δ G max Simulations Analytical O p t i m a l v o l t age ( m V ) Length of the shortest path (number of memristors) (b) (e) G ON G OFF (a) (c) (d)
FIG. 1. (a)-(c)-(d) Schematic of a memristor circuit at itssteady state for different applied voltages: 4 mV for panel(a), 1 mV for panel (c) and 10 mV for panel (d). The colorscorrespond to the conductance of each memristor. (b) Con-ductance versus time for the shortest path memristors (inred) and the longer paths memristors (in black). Note thatfor these ideal devices, the curves of all memristors in eachpath category superimpose. The difference ∆ G is indicatedby a double arrow. This corresponds to the graph depictedin panel (a), with an applied voltage of 4 mV. (e) Optimalvoltage (leading to the highest ∆ G ) versus the length of theshortest path, i.e., the number of memristors composing it.Each data point corresponds to a simulation on one of the1996 randomly generated square grid graphs. The solid linecorrespond to the analytical value NγτG ON . memristor connected to this node is measured. The nextnode on the shortest path is the one connected to thecurrent node through the highest conductance memris-tor. This process is repeated until the end node – thenode connected the ground – is reached. This methodallows the user to read the shortest path by probing onlya fraction of the memristors in the network, proportionalto the length of the shortest path. For this method toread the correct shortest path, it is required that at eachstep, the measured memristor belonging to the shortestpath has a higher conductance than the others. We thusdefine as metric of success, ∆ G , the smallest difference inconductance between a memristor on the shortest pathand a memristor outside the path but connected to thesame node of the path. The shortest path problem is suc-cessfully solved if ∆ G >
0. The highest possible successis ∆ G max = G ON − G OFF . In the successful example ofFigure 1(a), ∆ G (cid:39) ∆ G max , as observed in Figure 1(b).On the other hand, the results in Figures 1 (c) and (d)both exhibit ∆ G (cid:39) G is thehighest. Figure 1(e) shows that the optimal voltage isproportional to the length of the shortest path. In or-der to turn on a memristor, the time derivative of its x must be non-negative when x approaches 1, i.e, the volt-age across the device is greater than γτG ON . The lowestcontrol voltage to achieve this is NγτG ON , where N is thenumber of memristors on the shortest path. This opti-mal control voltage will turn on the shortest path butnot any longer paths. As observed in Figure 1(e), thismatches simulation results. Our simulations show thatusing a constant control voltage too far from the optimalcontrol voltage gives results like those in panels (c) and(d) Fig. 1, making such an appraoch impractical for ap-plications for which the length of the shortest path is notknown in advance. V. USING A VOLTAGE RAMP
We propose a method that does not require priorknowledge of the shortest path length. We leverage thefact that the shortest path is turned on (i.e., has the con-ductances of its memristors increase to G ON ) at lowercontrol voltages than longer paths, which suggest the useof a voltage ramp. Figure 2(a) shows the evolution of ∆ G with time as the voltage is increased. We observe a sharpincrease in ∆ G , which corresponds to the shortest pathturning on. The increase in conductance of the memris-tors on the shortest path creates an increase in the globalconductance of the circuit, shown in Figure 2(b). Thisincrease can be detected by measuring the current goingin and out of the circuit. The turning on of the short-est path corresponds to a sharp kink in the current, asshown in Figure 2(c), and thus a drop below zero in thesecond time derivative of the current, as shown in Figure2(d). When this drop is measured, the control voltage isturned off and the result of the computation is read outas described in Section IV. Note that here the evolutionof the system is stopped, contrary to the constant controlvoltage method where a steady state is reached.Figure 3 presents statistics on ∆ G obtained from sim-ulations and shows that this method, using a single volt-age ramp, can successfully find the shortest path in thou-sands of graphs of various sizes. In addition to the squaregrid based topology (blue bars), we have performed sim-ilar simulations on graphs with a small-world topology(red bars). Small-world networks stand in-between reg-ular and random graphs and describe many interesting Conductance (mS) Current ( m A) ( d )( c )( b ) d2Current/dt2 (nA/s2) T i m e ( s ) ( a ) D G / D Gmax
FIG. 2. Evolution with time of (a) the conductance difference∆ G/ ∆ G max , (b) the global conductance of the circuit, (c) thecurrent through the circuit and (d) the second time derivativeof this current. The voltage ramp starts from 0.1 mV andincreases at a rate of 0 . / s. The dashed line correspondsto the time of result detection. Count (cid:1)
G / (cid:1) G m a x S q u a r e g r i d S m a l l - w o r l d
FIG. 3. Histograms of the metric of success ∆ G/ ∆ G max forsimulations on randomly generated graphs, using the genericlinear model on 1996 square grid based graphs (blue bars)and 4797 small-world networks (red bars). The voltage rampstarts from 0.1 mV and increases at a rate of 0 . / s. problems and are discussed in more detail in Sec. VI.Here we have generated small-world networks of differentsizes and levels of randomness. For both topologies, allsimulations exhibit ∆ G well above zero, which validatesthis method.We test the validity of our method for realistic devices Count D G ( m S ) w i t h o u t v a r i a b i l i t y w i t h 1 0 % v a r i a b i l i t y
FIG. 4. Histograms of the metric of success ∆ G for sim-ulations on 658 randomly generated graphs following thesquare grid topology, using the experimentally obtained re-alistic model from Ref. 22 without variability (red bars) andwith 10 % variability on all device parameters (blue bars).With device variability, each simulation is done on one graphwith one set of parameters randomly chosen from a Gaussiandistribution around the nominal parameter value and of 10 %standard deviation. The set of graphs with and without vari-ability are the same. The voltage ramp starts from 0 V andincreases at a rate of 1 mV / s. by performing simulations on a square grid topology withthe realistic memristor model and parameters describedin Section III. Figure 4 shows that ∆ G is more widelyspread than for the generic linear model. This is dueto the fact that the realistic model produces smootherbehavior and a weaker reinforcement mechanism. How-ever, our approach remains valid, as all graphs exhibit∆ G > G than in the case with unique solutions. Qual-itatively, the system tries to turn all shortest paths onsimultaneously, which prevents the winner-take-all rein-forcement mechanism to take place properly. However,device variability makes one of the shortest paths intrin-sically more conductive and easier to turn on, which letsthe system choose this path over the others and turn it oncompletely, thus, increasing ∆ G . This is an interestingexample of how, in bio-inspired computing, features ofnanodevices usually seen as drawbacks can be beneficial. VI. SCALING OF THE TIME AND ENERGYCONSUMPTION
In order to evaluate the potential use of our method,we study how the time and energy required by the com-putation scale with the size of the graph. The energy wasestimated as the integral over time of the total currentthrough the circuit times the voltage across the circuit.Additional energy will be spent for detecting the currentsecond time derivative, reading the result and setting upthe circuit, but this is out of the scope of this study asdesigning the full architecture of the system would berequired.Figures 5(a) and 5(b) present the time and energy con-sumption versus the number of nodes plus the number ofedges for all our simulations with the generic model (cor-responding to Figure 3). This combination is a commonway to characterize the size of a graph. We observe nocorrelation. However, Figures 5(c) and 5(d) show thatthe time and energy consumption correlate strongly withthe length of the shortest path. Moreover, this scalingdoes not depend on the topology of the graph.These results indicate a key advantage of the pro-posed method. Conventional algorithmic methods typ-ically scale with the number of nodes and the number ofedges, because the different nodes and edges are exploredsequentially.
In the present hardware implementation,the current explores the entire circuit in parallel, whichmakes the time and energy consumption independent ofthe size of the graph. This method would be particu-larly efficient for large graphs with small shortest paths.Such graphs include small-world networks. Standing be-tween regularity and randomness, these are composed ofmany short range connections and a few long range con-nections. They exhibit high clustering and low shortestpaths. Small-world networks have been shown to de-scribe many systems with important applications, suchas power grids, the structure of the web, social media,and neural networks.
We investigate the effect of a realistic memristor modeland device variability by computing the time and energyconsumption corresponding to the simulations used inFigure 4. We observe that the scaling laws stay valid: asshown in Figure 6, the time and energy consumption de-pend on the length of the shortest path. The fact that thetime and energy consumption appear to scale indepen-dently of the topology or memristor model is promisingfor the ease of material implementation, as many typesof devices could be used, and for the breadth of applica-tions, as many types of graphs could be solved.
VII. OUTLOOK FOR A RECONFIGURABLEGRAPH SOLVER
In order for this approach to be of practical use, itmust be reconfigurable for different problems. Find-ing the shortest path between different nodes of the
Total time (s)
R e g u l a r g r i dS m a l l - w o r l d ( d )( c ) ( b )
Total energy (nJ)
N u m b e r o f n o d e s + n u m b e r o f e d g e s ( a )
S h o r t e s t p a t h l e n g t h
FIG. 5. Total time of the computation versus (a) the size(number of nodes plus number of number of edges) and (b)the length of the shortest path. Total energy consumed versus(c) the size and (b) the length of the shortest path. Bluecircles correspond to square grids and red squares correspondto small-world networks. Each symbol corresponds to one ofthe 1996 square grid and 4797 small-world graphs. ( b )
Time (s)
S h o r t e s t p a t h l e n g t h( a ) w i t h o u t v a r i a b i l i t y w i t h 1 0 % v a r i a b i l i t y
Energy (nJ)
S h o r t e s t p a t h l e n g t h
FIG. 6. (a) Total time of the computation versus the lengthof the shortest path. (b) Total energy consumed versus thelength of the shortest path. Each symbol corresponds to asimulation using the realistic memristor model either with-out variability (red full circles) or with 10 % variability oneach parameter (blue squares), on one of the 658 square gridgraphs. With device variability, each simulation is done onone graph with one set of parameters randomly chosen froma Gaussian distribution around the nominal parameter valueand of 10 % standard deviation. The set of graphs with andwithout variability are the same. same graph simply requires connecting the voltage sourceand ground to the new nodes. Modifying the graph isnon-trivial. Nearest neighbor connections could be im-plemented by complementary metal-oxide-semiconductor(CMOS) switches that are opened or closed to form thedesired graph, as proposed in Ref. 15. Longer range con-nections could be implemented by another set of mem-ristive devices, arranged in conventional crossbar topolo-gies. The hybrid CMOS/molecular (CMOL) architec-ture, proposed in Ref. 31, has been shown to give highdensities of connections, which could be useful here.
Building one physical system capable of implementingany graph problem is unrealistic. However, it would bepossible to have specialized chips for types of graph prob-lems. There would be one underlying hierarchical struc-ture appropriate to the graph type, as well as many re-configurable connections implementing specific problems.Such hierarchical architectures have been shown to be ef-ficient at simulating large artificial neural networks. VIII. CONCLUSION
We propose a scheme to find the shortest path in agraph problem using a circuit of memristors. It includesprocedures to detect and read the result. Our schemedoes not require prior knowledge about the shortest pathand has been validated on large number of graphs of var-ious sizes and topology. We have shown that this schemeworks for realistic device models and is robust to vari-ability. It scales with the length of the shortest path interms of time and energy consumption because of its in-trinsic parallelism. This is a key advantage compared toconventional algorithmic methods that scale with the sizeof the graph. In particular, it would be best for study-ing large graphs with small shortest paths, such as socialnetworks, neural networks, or power grids. These results are promising for hardware implementations of systemscapable of performing fast and energy efficient analysisof large graphs.More broadly, the field of swarm intelligence is rich,and implementing its concepts in hardware offers manypaths towards energy efficient computing. For example,it was shown that implementing a swarm intelligence al-gorithm of image edge detection with circuits of memris-tors consumes less energy than conventional methods. Exploring other swarm intelligence ideas and differentsubstrates to implement them is an exciting road towardslow energy cost systems that perform complex optimiza-tion tasks.
ACKNOWLEDGMENTS
The authors acknowledge A. Madhavan, M. Daniels,N. Zhitenev, J. McClelland, R. McMichael, S. Dushenkoand P. Shrestha for helpful comments and discussions.A.M. acknowledges support under the Cooperative Re-search Agreement between the University of Marylandand the National Institute of Standards and Technol-ogy, Center for Nanoscale Science and Technology, GrantNo. 70NANB10H193, through the University of Mary-land. This material is based upon work supported bythe National Institute of Standards and Technology Sum-mer Undergraduate Research Fellowship (SURF) Pro-gram under Grant No. 70NANB18H080. ∗ [email protected] E. W. Dijkstra, Numerische Mathematik , 269 (1959). R. K. Ahuja, K. Mehlhorn, J. Orlin, and R. E. Tarjan, J.ACM , 213 (1990). B. V. Cherkassky, A. V. Goldberg, and T. Radzik, Math-ematical Programming , 129 (1996). M. Thorup, J. ACM , 362 (1999). M. Thorup, J. Comput. Syst. Sci. , 330 (2004). M. L. Fredman and R. E. Tarjan, in (1984)pp. 338–346. R. Williams, in
Proceedings of the Forty-sixth Annual ACMSymposium on Theory of Computing , STOC ’14 (ACM,New York, NY, USA, 2014) pp. 664–673. E. Bonabeau, M. Dorigo, and G. Theraulaz, Nature ,39 (2000). M. Dorigo, M. Birattari, and T. Stutzle, IEEE Computa-tional Intelligence Magazine , 28 (2006). R. S. Parpinelli, H. S. Lopes, and A. A. Freitas, IEEETransactions on Evolutionary Computation , 321 (2002). C. Blum, Physics of Life Reviews , 353 (2005). M. N. A. Wahab, S. Nefti-Meziani, and A. Atyabi, PLOSONE , e0122827 (2015). F. Ducatelle, G. A. D. Caro, and L. M. Gambardella,Swarm Intelligence , 173 (2010). R. Poli, J. Kennedy, and T. Blackwell, Swarm Intelligence , 33 (2007). Y. V. Pershin and M. Di Ventra, Physical Review E , 046703 (2011). Y. V. Pershin and M. Di Ventra, Physical Review E ,013305 (2013). Y. V. Pershin and M. D. Ventra, Neural Processing Letters , 265 (2016). Z. Pajouhi and K. Roy, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems , 1774(2018). D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S.Williams, Nature , 80 (2008). M. D. Pickett, D. B. Strukov, J. L. Borghetti, J. J. Yang,G. S. Snider, D. R. Stewart, and R. S. Williams, Journalof Applied Physics , 074508 (2009). S. Menzel, M. Waters, A. Marchewka, U. B¨ottger,R. Dittmann, and R. Waser, Advanced Functional Ma-terials , 4487 (2011). T. Chang, S.-H. Jo, K.-H. Kim, P. Sheridan, S. Gaba, andW. Lu, Applied Physics A , 857 (2011). M. H. R. Lankhorst, B. W. S. M. M. Ketelaars, andR. a. M. Wolters, Nature Materials , 347 (2005). D. Ielmini and A. L. Lacaita, Materials Today , 600(2011). S. Lai, in
IEEE International Electron Devices Meeting2003 (2003) pp. 10.1.1–10.1.4. R. Bez, in (2009) pp. 1–4. D. J. Watts and S. H. Strogatz, Nature , 440 (1998). R. Albert, H. Jeong, and A.-L. Barab´asi, Nature , 130 (1999). D. S. Bassett and E. Bullmore, The Neuroscientist , 512(2006). G. Palla, I. Der´enyi, I. Farkas, and T. Vicsek, Nature ,814 (2005). D. B. Strukov and K. K. Likharev, Nanotechnology ,888 (2005). D. Strukov and A. Mishchenko, in
Proceedings of the Con- ference on Design, Automation and Test in Europe , DATE’10 (European Design and Automation Association, 3001Leuven, Belgium, Belgium, 2010) pp. 661–666. A. Madhavan, T. Sherwood, and D. B. Strukov, IEEETransactions on Very Large Scale Integration (VLSI) Sys-tems , 1 (2018). R. M. Wang, C. S. Thakur, and A. van Schaik, Frontiersin Neuroscience12