Asymptotic Behavior of Memristive Circuits

The interest in memristors has risen due to their possible application both as memory units and as computational devices in combination with CMOS. This is in part due to their nonlinear dynamics, and a strong dependence on the circuit topology. We provide evidence that also purely memristive circuits can be employed for computational purposes. In the present paper we show that a polynomial Lyapunov function in the memory parameters exists for the case of DC controlled memristors. Such a Lyapunov function can be asymptotically approximated with binary variables, and mapped to quadratic combinatorial optimization problems. This also shows a direct parallel between memristive circuits and the Hopfield-Little model. In the case of Erdos-Renyi random circuits, we show numerically that the distribution of the matrix elements of the projectors can be roughly approximated with a Gaussian distribution, and that it scales with the inverse square root of the number of elements. This provides an approximated but direct connection with the physics of disordered system and, in particular, of mean field spin glasses. Using this and the fact that the interaction is controlled by a projector operator on the loop space of the circuit. We estimate the number of stationary points of the approximate Lyapunov function and provide a scaling formula as an upper bound in terms of the circuit topology only.


Introduction
The study of memristors has become recently an area of interest [1][2][3] for a variety of reasons. This is not only due to the fact that memristors are considered by many the fourth circuit element, but also due to their potential applications and their interesting dynamics. In particular, it has been understood from a purely computational theory perspective that the combination of memristive and CMOS components leads to universal computing machines, called memcomputers [4]. From a theoretical perspective, circuits made of memristors have been shown to exhibit non-trivial dynamics which can in principle serve for computational purposes [6][7][8][9]. Memristors are passive components which can be thought of as a time varying resistance sensitive to the either the current or the voltage, which in turn depends on a dynamical internal state variable. The main characteristic of a memristor is a pinched hysteresis loop in the Current-Voltage diagram when controlled in alternate voltage [11][12][13]. The aim of this paper is characterizing the asymptotic behavior of purely memristive circuits (i.e. only memristors) for general circuit topology. An open question which has not insofar been answered at a theoretical level is what is the role of the circuit topology in the relaxation process of memristive circuits. It is in fact thought that the topology plays a role in the optimization capabilities of these circuits with memory (also called information overhead by Di Ventra and Traversa [4]). We provide a precise answer for the simpler class

Memristive circuits
Before we introduce the dynamics for a generic circuit, we first briefly discuss the type of memristors under scrutiny [11]. Specifically, we consider memristors whose internal dynamics (the parameter w) depends on the current only, and in which the resistance depends linearly on an internal parameter w, and satisfy the linear relation R(w) = R on (1 − w) + R o f f w, with 0 ≤ w(t) ≤ 1. We consider in particular the time evolution of a single memristor with diffusive dynamics [14,15] and for an applied voltage S: which shows that a competition between drift and decay occurs. This type of memristors has been considered for machine learning applications. With our convention, R on and R o f f are the limiting resistances for w = 0 and w = 1 respectively (R o f f > R on > 0), and α and β are constants which set the timescales for the relaxation and excitation of the memristor respectively. 1 In a recent paper [16] it has been observed that mean field theory and techniques from statistical physics can be applied to study a specific circuit topology. In that paper, it was also noted that a Lyapunov function exists, and that zero temperature mean field theory provides a good estimate for the average asymptotic dynamics for a single mesh of memristors. In the present paper we extend some of these results to a more general class of memristive circuits, which sheds new light onto this type of nonlinear systems. Specifically, in [6] the following differential equation was derived for a generic but purely memristive circuit: 1 It is worth mentioning that our model is different from the one introduced in [11] as it has not only opposite polarity, but also a decay function which is slightly different. where S is the vector of voltage sources in series to each memristor, while I is the identity matrix. It is immediate to observe that the nonlinearity is controlled by the parameter ξ = . It is important to note that Ω and W are matrices, while I is the identity matrix. W is the diagonal matrix with the memristor memory values w i (t) on the diagonal, meanwhile Ω = A t (AA t ) −1 A is a projector operator on the fundamental loops of the circuit, being A the cycle matrix of the directed graph representing the circuit flows 2 . The richness of the dynamicsl of this equation has been characterized in [6], while the locality properties of Ω in [17]. The important observation that we anticipate is that the dynamical eqn. (2) has at least one Lyapunov function which we have derived analytically. While the derivation is provided in the Appendix, here we provide only its functional form: from which a few key facts can be immediately observed. For instance, the role of Ω is the one of a coupling matrix, as one would expect, and that it is not quadratic in the internal memory variables. Eqn. 3 satisfies d dt L(W) ≤ 0 whenever M i = −2αξ ∑ j =i W j Ω ji W i is small (thus for very weakly interacting memristors, as it depends only on the offdiagonal terms) and d dt L(W) = 0 ↔ d dt W = 0 (see the details in the Appendix). One can obtain a more precise bound on the derivative of the Lyapunov function. Let us define ||Ω S|| 2 = N 2 s 2 (N). We have proved in the Appendix that if then d d L(W) < 0, whereΩ = max ij |Ω ij |. It is interesting to note that the order parameter s(N) αβ is a generalization of the one found in [16] and which controlled the asymptotic state of the effective mean field circuit. Here again it plays a role. Note that ξ 2 (1 + ξ)Ω 2 is what controls the height of the parabola, and intuitively the larger then nonlinearity the higher the voltage has to be for the system to relax according to the Lyapunov function. Thus, eqn. (4) establishes a dynamical phase diagram.
We stress that the dynamics is constrained in the hypercube [0, 1] N , where N is the number of memristors, and thus the Lyapunov function above when a memristor reaches the boundaries of this domain. However, if we now use the fact that the dynamics is controlled with constant voltage each memristor will eventually reach the asymptotic value 1 or 0. In this asymptotic limit, we note that the 2 The matrix A has only −1, +1, 0 and zero values. Given an orientation of each loop and edge, then A iβ has dimensionality L × N, where N is the number of memristors and L the number of fundamental loops. If a memristor β has the same orientation of a loop i, then A βi = 1, and −1 if opposite. A βi = 0 if the memristor β does not belong to the loop i. The matrix Ω is not diagonal because edges of a graph can belong to more than one loop.  function above can be interpreted as a spin-like model asymptotically. 3 In this case, the Lyapunov function reads where we have used the fact that W n i = W i if W i = 1, 0. An effective external field h i = α 2 + αξ 3 Ω ii − 1 β Ω ij S j emerges in this asymptotic approximation. We also test numerically both the Lyapunov function of eqn. (A1) and eqn. (5) in Fig. 1. Meanwhile in Fig. 1 (a) we plot the evolution of each memory element, in Fig. 1 (b) we plot the evolution of both the Lyapunov function and the asymptotic approximation of it, which remains close to the exact one. On the other hand, in Fig. 2 we show the difference between the obtained Lyapunov function and the asymptotic one as a function of time. This mapping is reminiscent of the case of continuous neuronal networks introduced by Little [19] and then Hopfield in a series of important papers [20,21]. The Little-Hopfield model has sparked interest from the Statistical Physics community since the very beginning [22]. In the past years this particular line of study has been subject of scrutiny by many experimental groups [23][24][25]. The key difference is that in each case, these were studied in conjunction with ordinary and active (CMOS) electronic components to build Hopfield learning networks. We argue instead that memristive circuits per se form a special kind of Hopfield network defined by the Lyapunov function above, without the need of extra components.

Functional complexity for random circuits
We now study random circuits, in particular for the function complexity of eqn. (5). First, this will allow us to show an interesting connection to the field of disordered systems, and in this approximation provide estimates of the complexity of the Lyapunov function. In fact, minimizing quadratic function with discrete variables is in general a hard problem to solve. We are however in the situation in which the problem at hand is in nature both continuous in time and in the variables (being these between 0 and 1), but naturally provides an answer for the more complicated case usually associated to finding the ground state of a frustrated spin systems. Despite the fact that we are not able in the present paper to provide a complete answer to which optimization class a system of purely memristive circuits belongs to, we try to provide some answers to the questions above using some techniques introduced to study random polynomials [26,27]. This will be important in light of the fact that we have shown that memristors perform a local optimization, rather than a global one, and that the circuit constraints should enter somehow. 3 In [16] this phenomenon has been partly explained by the emergence of unstable fixed points in the dynamics (although only in the mean field approximation). For the specific case of the eqn. (2) however, it is necessary to go beyond this approximation, for which each memristors will have (intuitively) a different fixed point associated to it. This phenomenon, which is circuit dependent, is associated to the structure of the couplings Ω ij , and will be studied elsewhere.
The number of stationary points of a generic multi-varied function L(W) can be estimated by: Thus, if we aim to consider the expected number of stationary points #, then it does make sense to consider the (quenched) average of the quantity in eqn. (6) for the distribution for · P(Ω) with respect to the coupling matrix Ω. One important observation is that because Ω is a projector, det(∂ W i ∂ W j L(W)) can be calculated exactly without knowledge of the elements distribution, but only using the fact that Ω is a projector (details in the Appendix). We find that the average number of stationary points can be split in the form: where L is the number of fundamental loops of the circuit, and where: depends on the distribution of the elements of the matrix Ω.
We first need, then, to understand the distribution P(Ω) for random circuits. Since the analytical calculation of such distribution goes beyond the scope of this paper, we perform such analysis numerically. We first generate random circuits with Erdos-Renyi graph above the percolation threshold (p = 0.7). We are interested, in particular, in the scaling with respect to the number of memristors N. For fixed Ω, the distribution P(Ω) of off-diagonal elements is shown in Fig. 3 (a). We observe that although this seems to be unimodal, a careful look shows that it is not. We observe in Fig. 3 (b) that there is a non-zero probability of having elements not distributed around zero, but that these are orders of magnitudes smaller than the central distribution. While it is not surprising that the distribution is not completely random (thus, there are correlations among elements) since these are constrained by Ω 2 = Ω, it is surprising to note that we can roughly approximate this distribution with a unimodal one. We are in particular interested in how the width of the distribution scales with the number of memristors. This can be calculated exactly if we know the scaling of the diagonal elements, as we will see below. In Fig. 4 we show how the diagonal elements of Ω scale with N. We fit first the diagonal elements, precisely 1 − Ω ii , showing that it scales as a power law to a good approximation. If we perform a fit using the functional form Ω ii ≈ 1 − c N α , we obtain the best fit values c = 1.74557 ≈ √ 3 and find that α = 0.5043 ≈ 1/2. We will, since now on, consider these values in what follows. It is important to note that the scaling above is sufficient to obtain a scaling for the off-diagonal elements as well. We define the matrix of off-diagonal element G, Since we have that Ω 2 = Ω, it is easy to see that G must satisfy the equation as well. We can at this point solve for an eigenvalue equation for G, and obtain λ = { 3 N , −1 + 3 N }. Since this must be a projector for large values of N, we must have Ω ≈ I + 3 N Q, where G = 3 N Q, with Q again a projector. Thus, the scaling approximation obtained for random Erdos-Renyi type circuits, we This is the scaling we will use in the following. While in what follows we show only the intermediate results, all the exact calculations are provided in the Appendix. We can in fact a this point calculate Z in eqn. (6) for random circuits. We assume that where σ is an effective width for N = 1. In this case, calculating Z can be done by means of Gaussian integrals. In this case, we obtain, in the limit αξ 1: Figure 5. Obtained minimum energy using the minimum over 100 random binary values (blue crosses), roghly 8000 steps of exponential annealing (green pluses) and 60 time steps using the memristor dynamics (red circles). We generated each sample with approximately 800 − 900 memristors, fixing the percolation parameter to p = 0.7, with ξ = 10, α = 0.1 and β = 1. The numerical integration was performed using a simple Euler-Newton integration with integration step dt = 0.1, and total number of steps T = 1000. We truncated the Metropolis annealing to 10 * N time steps, with N the number of memristors, and with an annealing rate λ = 0.995, and an exponential annealing law for the temperature given by Temp(k) = 100 * λ k .
We now use Euler's relation L = N − V + χ, where V is the number of vertices, χ the Euler characteristic and use the Erdos-Renyi relation where p(V) is the edge probability. Now we have that, since χ > 0 and only topological, while L scales linearly for large N in the number of memristors. We now observe that given ρ = √ 3 √ πσ , we have that in the limit N → ∞, # goes to infinity for ρ > 1, while it goes to zero for ρ ≤ 1, as lim We thus obtain the result that for σ > √ 3 √ π the Lyapunov function seems to have a large number of stationary points, and thus a hint of the hardness of the minimization problem.
It is easy at this point to provide an approximate mapping between the asymptotic Lyapunov function and the Hamiltonian of a mean-field spin glass of the Sherrington-Kirkpatrick type with an external field: where we used the mapping and disregarded an unimportant constant which only shifts the function. The result above is a generalization of what has been found in [18] in the limit ξ 1, e.g. the fact that the dynamics of memristors follows a constrained gradient descent (Rosen projections). Moreover, it provides a physical system to test experimentally the physics of mean field spin glasses [28].

Asymptotic state recollection and combinatorial optimization
One important question for memristive circuits is what is the asymptotic state they reach. As we have argued, this can be answered by looking at the minima of the Lyapunov function for the case of constant voltages. Although this is a hard problem to solve in its full generality, we can use some approximate analytical formulae to provide an answer in at least for some region of the parameters. First, we consider the case α = 0. It is easy to see that we can integrate the equation [18] to obtain: where c is an integration constant. Eqn. (12)  . For ξ > 0 and in the asymptotic limit, lim t→∞ W 2 ≈ W. Using this approximation, we obtain: We then can obtain the asymptotic formula: where we note that in the last equation we have used the fact that 1 + ξ/2 > 0 if 0 < R on < R o f f . However, we note that the exactness of the asymptotic behavior depends only on ξ and not on α. The results on how well eqn. (13) predicts the asymptotic behavior of the circuit as a function of ξ are shown in Fig. 6. Despite our heuristic method and approximation, the obtained values are a good approximation for small values of ξ of the real system. Using the formula above, we now provide the connection with the state recollection of neural networks. We note that the projection operator Ω can be written as [17] Ω ij = ∑ L l=1Ã l iÃ l j , wherẽ A l i is the orthonormalized loop matrix of the circuit. From the theory of neural networks we are aware that the number of stored patterns equals the number of independent eigenvectors of the interaction matrix Ω, which is purely topological, as argued in [18]. In the case of a purely memristive circuit the number of independent memory units number is constrained by the topology of the circuit, and it depends on the number of fundamental loops [18]. As we have noted before, the asymptotic behavior of the circuit can be approximated by (where we reinstate α) where we realize that the source vector S can act as the recollection mechanism. If S = ∑ l ρ lÃ l , for a certain proportionality constant ρ, then (Ω S) i = ∑ l ρ lÃ l i . Thus, the asymptotic configuration of the memristors will be given by the approximate value: which shows that stored patterns are indeed in the loop basis. The analysis above is reminiscent of the asymptotic state recollection in the Hopfield-Little model. We now ask the converse question of how close are the asymptotic states to minima obtained from the Lyapunov functional. As we have seen in the previous sections memristors can be used to provide solutions (possibly sub-optimal) to hard combinatorial problems. This claim should be taken, if not with skepticism, with some care. First we note that L(W) is not positive definite (although it is bounded), and thus we have an asymptotically local stable equilibrium rather than a global one. We are interested in providing some benchmarks on how well memristors find the minimum of the function L a (W) compared to other optimization techniques. Given 100 samples of random circuits based on Erdos-Renyi random graphs with p = 0.7, in Fig. 5 we show the energy attained using memristors in 60 time steps (red dots), using a Metropolis algorithm with exponential annealing (blue) with over 8000 steps, and using the minimum of 100 random configurations as a reference (blue dots). We see that despite the fact that memristors do not seem to reach the absolute minimum, the energy is closer to the Metropolis result than to the random configuration one. Although it is easy to see that the dynamical equations do not reach the Metropolis attained minimum, it did take less than 100 iterations for the memristive system to converge using a simple Euler-Newton first order integration. This already shows that, as anticipated, memristors perform a local optimization in this regime. In the next section, however, we apply the algorithm to a specific dataset, where this optimization has some advantage over the exponential annealing.

An application: minimization of a function
We have argued that memristors' dynamics can be thought, when these are connected in a circuit, as reaching an asymptotic state which is among the minima of a certain Lyapunov function. Albeit real memristors can have a dynamics which is far from ideal, we would like to consider if these could in principle be used for any sort of purpose which goes beyond electronics. In fact, note that (2) can be used as a quick meta-heuristic method to search for the minimum of a function, or to be given as input to other more complicated and efficient optimization techniques. In realistic circuits, the quadratic form matrix Ω which occurs in the Lyapunov function is a projector on the cycle space of the graph. Thus, problems for which analog memristive circuits can be used are naturally those that can be embedded in the cycle basis of a graph [34]. However, in principle we can simulate the system also for matrices Ω that are not necessarily projectors, as we have shown and as we discuss further below. We can thus perform a preliminary test of the applicability of analog memristive circuits in a problem of optimization, and see how the circuit performs compared (for instance) to simulated annealing.
We thus consider a simple application of practical importance, e.g. studying the problem of investment in a set of assets. For this purpose, we use the 225 Nikkei dataset which is used for benchmarking heuristic optimization algorithms [35], available at [36]. We can use in fact a memristor-inspired optimization scheme, taking advantage of the fact that it works beyond projector operators Ω. We mention that the proof of the Lyapunov function relies only on the matrix I + ξ ΩW+WΩ 2 being positive definite, and thus applies also to any positive matrix and to some non-positive ones (see SI for details) for ξ small enough. Although we cannot use these in a hardware circuit, we can solve the differential equation numerically to infer a minimum heuristically. Using this idea, we ask in which of the assets of portfolio we should invest in (a yes-no answer). The dataset is composed of 225 assets, including returns and the covariance matrix, and can thus use the combinatorial Markowitz functional [37]. If we use the fact that maximizing a quadratic function is equivalent to minimizing the function with a minus sign in front, we can write an equality between the Markowitz function for binary variables and the memristive one: where p is a trade-off parameter, r are the returns and Σ is the covariance matrix. From imposing the equality between the two functionals, we observe the necessity of imposing Ω = Σ, and thus We can thus obtain the source vector to try to force the system towards the right minimum: where (η) i = Σ ii , e.g. the variance of each return. We have the freedom to choose arbitrarily β, but either ξ or α would be fixed. We observe that we need to invert the covariance matrix, which is a slow but Figure 7. Maximization of the return for the Nikkei 225 dataset using a Monte Carlo procedure. We compare the results between the best of 100 Simulated Annealing procedures (the green line in the plot ≈ 9.15) versus with the memristive optimization proposed in this paper (black lines), with identical initial conditions. The simulated annealing was conducted with an initial temperature T = 100, and an exponential annealing rate λ = 0.999 (and effective time steps N=500). We have then used the final states of the Memristive Optimization procedure as an input to a simulated annealing procedure with lower temperature T 0 = 0.025 and identical annealing rate, which obtained the absolute maximum of the return.
polynomial in the number of variables. Moreover, we see that the harder the problem to solve (the norm of the inverse of Σ), the smaller β has to be chosen, from which we infer the slowness of the process.
We are now in the position to test the memristive minimizationi against an exponential annealing process, which is known to perform poorly in the combinatorial setting. The results for the case of the Nikkei dataset are shown in Fig. 7, comparing the Metropolis annealing to the case of the computational time and in the final value obtained. We performed three different tests. First, we randomized the initial states, and performed 100 Simulated Annealing procedures [38], with initial temperature T = 100, an exponential annealing rate λ = 0.995 (T k = λ k T 0 ) and a number of time steps of N T = 500 * N, where N was the number of assets 100. On the Nikkei 225 dataset, the maximum obtained with the simulated annealing was r max ≈ 9.15. With the identical initial states, the memristive optimization has obtained a maximum of r max ≈ 14.23. We have then used the final state of the memristive optimization output as an input to an annealing procedure, but with a lower initial temperature T 0 = 0.025, thus effectively fast-forwarding the annealing procedure. This combined optimization obtained the absolute maximum, albeit of only less than a percentage better than the memristive optimization, of r max ≈ 14.40.

Concluding remarks
In the present paper we have analyzed the asymptotic dynamics of memristive circuits, showing that an approximate Lyapunov function for the dynamics exists for memristors with slow decay. Because of the properties of memristors, we have argued that asymptotically the Lyapunov function of a memristive circuit can be written as a quadratic function on spin-like variables, and have connected the coupling matrix to the projector operator on the loop space of the graph [10]. This result shows a direct connection between purely memristive circuits and Hopfield networks, which we argue minimize a similar type of Lyapunov function. The internally stored patterns are in this case the cycles of the graph. This is interesting for various reasons. Insofar it had been argued that a certain degree of external control was necessary in order to perform computation or use memristors as a memory device, either via the use of CMOS or the introduction of capacitors into the network. This paper provides sufficient evidence for the use of memristors in their completely analog regime for computational purposes. This comes at a cost, which is the embedding of the problem of interest in the cycle basis of the graph, which can be a rather nontrivial problem when it is solvable.
The Lyapunov function for the dynamics of memristive circuits can then be used for various purposes. We have first focused on the complexity of these functions in the case of random circuits. We have provided numerical evidence and argued that the couplings scale as in the case of mean-field spin glasses, i.e. the Sherrington-Kirkpatrick model, although only if neglecting correlations between matrix elements of lowest order. Using this approximation we have provided approximate formulae for the number of stationary points, provided some topological considerations. This result explains the importance of the initial conditions for the asymptotic state of the memristors. In fact, it is known (for instance in the case of a glassy system) that for a system which has a large set of local minima there is a large sensitivity to the initial conditions.
As a test-bed of these ideas, we have used the memristor dynamics to see how well (or bad) the dynamics leads to the minimization of a certain quadratic functional We have compared the minimization of a combinatorial problem using both simulated annealing and the memristive dynamics introduced in this paper. Although the memristor dynamics could not reach a global minimum, the speed of convergence and the closeness to the Metropolis result suggest the use of these mixed dynamical-combinatorial algorithm to provide quick answers to combinatorial problems on spin-like variables (or alternative 0-1 variables, thus a quadratic unconstrained binary optimization). We do not claim these answers to be optimal, but we find nonetheless interesting that a rather simple procedure like the one we did performed remarkably well compared to simulated annealing.
The main result of this paper however remains the Lyapunov functional for the study of the dynamics of memristors, and to understand the relaxation properties of these circuits (at least when these are controlled with constant voltage). We have shown a direct connection between optimization and memristive circuits, as advocated by other authors, for instance [30][31][32]. It is inspiring to think that the minimization of the functional depends on the balance between reinforcement-decay properties of the memristive dynamics, along the lines of other similar algorithms [33] which use collective dynamics as heuristic optimization methods.
Thus, this paper provides needed background work to understand the dynamics of circuits with memory, their sensitivity to initial conditions and their use for computational purposes in the fully analog regime.

Appendix A Appendix: Properties of L(W). Appendix A.1 Derivation of L(W)
The Lyapunov functional provided in the paper is given by: In order to derive the functional above, we use the differential equation for the memristive dynamics: for the case of constant voltages S. At this point we can obtain the functional form of the derivative of the Lyapunov function: Where we assumed that ∂ t L(W) = 0, as S is constant in time. We also introduced the quantity M i = −2αξ ∑ j =i W j Ω ji W i . The variation of the functional is provided in the subsection below in detail. First we prove that I + ξΩW is a positive definite matrix. We use the matrix similarity for eigenvalues, Since Ω is a projector and is semipositive, then we have that for any non-zero vector b, b t (I + ξ √ WΩ √ W) b > 0. We now observe that the first term in eqn. (A3) is always negative, meanwhile the second term is always positive, since the inverse of a positive operator is always positive. However, meanwhile the first term does not scale with any parameter, we have that the second term scales roughly as α 2 ξ 2 . Also, the min norm of the operator I + ξΩW is of order one, which implies that the sup norm of its inverse is also of order one, for arbitrary values of ξ. As a comment, we see that up to this point, this is enough to prove that if Ω is diagonal the second term disappears, as M = 0. This completes the proof that d dt L( W) ≤ 0 for the case of Ω diagonal, which is trivial. Also, note that we have d dt L( W) = 0 if and only if d W dt = 0. This confirms the fact that L( W) is a Lyapunov function in this regime, and thus applies also for weakly interacting memristors.
We are interested in the more general case of arbitrary Ω.
whereΩ = sup|Ω ij | is the maximum absolute value of the elements of Ω ij . On the other hand, we have and thus We thus find that if the Lyapunov function we introduced has negative derivative always. Let us now define ||Ω S|| 2 = N 2 s 2 (N). We can rewrite the inequality as Because s(N) ∼ O(1) in the physical case (e.g. the single voltage element does not scale with the size of the system), the bound above gives a meaningful relationship between the nonlinearity parameter ξ and the mean field dynamical properties of the system, s αβ , which had been already introduced in [16] in a mean field toy model of memristive dynamics. It is interesting to note that if we define s αβ = q, the equation above forms a parabola and an inequality of the form which characterizes the area below the curve. Thus, via the bound we see that below the parabola the system is going into a minimum of the Lyapunov function. This shows what we had anticipated, e.g. the fact that the Lyapunov function we defined has a negative derivative, and because we have that d dt L(W) = 0 ↔ d dt W = 0 we have the required property for ξ small enough. Also, being L(W) defined on a compact and being smooth, we know it must be bounded from below, which concludes the proof. Clearly, how proof has a main fallacy: it applies to continuous dynamical systems with continuous derivatives. If certain memristors reach the boundary values, discontinuities appear in the dynamics. Simulations however show that the Lyapunov functional we obtained is an upper bound to the one obtained numerically and including the constraints.

Appendix A.2 Variation
Let us now calculate the variation of the functional of eqn. (A3). Again, we consider the functional: We have that: where on the last line we used the equations of motion, and defined M i = −2αξ ∑ j(i =j) W j Ω ji W i . For random networks, This mapping is reminiscent of the case of continuous neuronal networks introduced by Little [19] and then Hopfield in a series of important papers [20,21]. The Little-Hopfield model has sparked interest from the Statistical Physics community since the very beginning [22]. In the past years this particular line of study has been subject of scrutiny by many experimental groups [23][24][25].
The key difference is that in each case, these were studied in conjunction with ordinary and active (CMOS) electronic components to build Hopfield learning networks. We argue instead that memristive circuits per se form a special kind of Hopfield network defined by the Lyapunov function above, without the need of extra components.

Appendix B Complexity of the Lyapunov functional via Kac-Rice formula
This section provides the details for the average number of stationary points of the Lyapunov functional above in the asymptotic regime. Let us now consider the following average: where L( w) is the asymptotic functional in eqn. (A3). We can see that and thus We are interested in an asymptotic upper bound on the number of local minima of the function L(w). We first proceed by calculating the determinant. If N is the number of memristors, then det(αξ Ω ij − Ω ii δ ij ) = (αξ) N det Ω ij − Ω ii δ ij . In the scaling observed above, we have noted Using this scaling, we have: det if A is invertible, and write: is also a projector on the subspace, we have: It is now easy to see that since v l · v l = 0 if l = l , we have that and thus we have the remarkable result that in the limit described above, we have which is independent from the distribution of Ω's elements, but only on the scaling we used. We thus obtain a rough scaling of the number of minima, given by: Up to here no assumption has been made on the distribution of Ω ij , if not the scaling observed numerically. In order to evaluate Z( S, Ω), however, assumptions have to be made. We now focus on evaluating the term: which can be written as: We will now introduce again the scaling Ω ii = 1 − 1 √ N , and Ω ij =  ) and b i (w) = αξw + 1 β S i , we can write the integral as: Z(S, P(Ω)) = 1 2π We now rescale σ 2 → σ 2 /N in order for σ to be physical in the limit N → ∞. We thus obtain, in the approximation in which S i = S, Z(S) becomes: (1) αξ . Thus, we have the final formula: which shows two regimes. If σ < √ 3/π ≈ 0.977205, the number of stationary points increases exponentially.