Critical branching processes in digital memcomputing machines

TL;DR

This paper demonstrates that digital memcomputing machines naturally self-organize into a critical state characterized by avalanche distributions following a Borel distribution, which enhances their efficiency in solving complex combinatorial problems.

Contribution

It introduces a mean-field theoretical prediction of avalanche size distribution in memcomputing machines and confirms it through numerical experiments on random 3-SAT instances.

Findings

Avalanche sizes follow a power-law distribution with an exponent around 1.5.

Memcomputing machines self-tune to a critical state during problem solving.

Numerical results align with the theoretical critical branching process model.

Abstract

Memcomputing is a novel computing paradigm that employs time non-locality (memory) to solve combinatorial optimization problems. It can be realized in practice by means of non-linear dynamical systems whose point attractors represent the solutions of the original problem. It has been previously shown that during the solution search digital memcomputing machines go through a transient phase of avalanches (instantons) that promote dynamical long-range order. By employing mean-field arguments we predict that the distribution of the avalanche sizes follows a Borel distribution typical of critical branching processes with exponent $\tau= 3/2$. We corroborate this analysis by solving various random 3-SAT instances of the Boolean satisfiability problem. The numerical results indicate a power-law distribution with exponent $\tau = 1.51 \pm 0.02$, in very good agreement with the mean-field analysis. This indicates that memcomputing machines self-tune to a critical state in which avalanches are characterized by a branching process, and that this state persists across the majority of their evolution.

Equations

$$C_n(v_i,v_j,v_k)=\frac{1}{2}\text{min}(1\pm v_i,1\pm v_j,1\pm v_k),$$

$$\begin{array}{r}\frac{d}{dt}{v}_i=-\sum _n{x}_{s,n}{x}_{f,n}{G}_{n,i}(v_i,v_j,v_k)+\\ (1-x_{f,n})R_{n,i}(v_i,v_j,v_k),\end{array}$$

$$\frac{d}{dt}{x}_{f,n}=\beta (x_{f,n}(1-x_{f,n}){)}^{1/2}(C_n(v_i,v_j,v_k)-\gamma ),$$

$$G_{n,i}(v_i,v_j,v_k)=\frac{\partial }{\partial v_i}(1\pm v_i)\text{min}[(1\pm v_j),(1\pm v_k)],$$

$$\begin{array}{rl}{R}_{n,i}(v_i,v_j,v_k)& =\\ & \{\begin{array}{ll}{v}_i\pm 1,& C_n(v_i,v_j,v_k)=\frac{1}{2}(1\pm v_i)\\ 0,& \text{otherwise}\end{array}\end{array}$$

$$\frac{d}{dt}{x}_{s,n}=\alpha x_{s,n}(x_{f,n}-C_n(v_i,v_j,v_k)),$$

Full text

I Introduction

Unconventional computing paradigms that employ physical properties to compute specific problems are emerging as an important research direction in Physics Nielsen and Chuang (2010); Kozma et al. (2012); Nawrocki et al. (2016). One such paradigm is memcomputing Di Ventra and Pershin (2013); Di Ventra and Traversa (2018), which employs time non-locality (memory) to both process and store information on the same physical location. The digital version of this paradigm has been introduced to specifically tackle combinatorial optimization problems Traversa and Di Ventra (2017). Digital memcomputing machines (DMMs) can be physically realized as non-linear dynamical systems whose point attractors represent the solutions of the problem to be solved.

Since DMMs are non-quantum systems, their equations of motion can be efficiently integrated numerically. Results from these simulations have already demonstrated that DMMs perform orders of magnitude faster than traditional algorithmic approaches on a wide variety of combinatorial optimization problems Traversa et al. (2018); Di Ventra and Traversa (2018); Manukian et al. (2018); Traversa and Di Ventra (2018); Sheldon et al. (2018).

Subsequently, by employing topological field theory Di Ventra et al. (2017), it was shown that the physical reason behind this efficiency rests on the dynamical long-range order that develops during the transient dynamics where avalanches (instantons in the field theory language) of different sizes are generated until the system reaches an attractor Sheldon et al. (2018). The transient phase of the solution search of DMMs therefore resembles that of several phenomena in Nature, such as earthquakes Bak and Tang , solar flares Lu and Hamilton (1991), quenches Pruessner (2012), etc. Since all these phenomena show scale-free properties in the probability distribution of the avalanche sizes, it is natural to ask whether DMMs would also share this property. In this paper, we indeed show that the transient dynamics of a DMM are characterized by a critical branching process. We first provide a general mean-field analysis to argue that the probability distribution of the avalanche sizes should be a critical Borel distribution with exponent $\tau=3/2$ Zapperi et al. (1995), irrespective of the problem to solve, namely it is an intrinsic feature of DMMs. We then support these results with numerical simulations of DMMs’ equations of motion applied to the solution of Boolean satisfiability (SAT) instances. We have chosen to work with randomly-generated, satisfiable 3-SAT benchmark instances precisely to ensure that any feature produced by our analysis is a feature of the dynamics of DMMs, rather than a feature of the SAT instances solved.

Random 3-SAT belongs to the class of propositional logic in which a formula of Boolean variables must hold true for the problem to be satisfiable Bradley and Manna (2010). Propositional variables appear as literals in the formula, where a literal is a variable or its negation. Satisfiability problems are traditionally represented in “conjunctive normal form” (CNF), i.e., a conjunction (AND) of disjunctions (OR) of literals Arora and Barak (2009). A disjunction of literals is referred to as a clause. Therefore, a 3-SAT problem is one in which all clauses contain three distinct literals, of which none is a negation of the others.

A CNF formula has a simple Boolean circuit representation Arora and Barak (2009). An example for a 3-SAT with three clauses is reported in Fig. 1. A DMM that solves the 3-SAT, say, of Fig. 1 can then be realized as an electrical circuit with memory (see Eqs. (2), (3), and (6) below) where each variable of the 3-SAT problem is represented by a voltage (we represent with $+1$ the logical 1 and with $-1$ the logical 0, in arbitrary units), and each traditional OR gate of Fig. 1 is replaced by a “self-organizing” OR gate Traversa and Di Ventra (2017), namely one that always attempts to dynamically satisfy the logical OR truth table at its terminals. Since the problems we are seeking to solve are satisfiable, the “out” terminals of the CNF formula in Fig. 1 are set to be logically true, hence the voltages at those terminals are kept fixed at +1.

II Mean-field analysis

With these preliminaries we can now discuss the transient dynamics of DMMs and argue that the size of the generated avalanches (of the voltages at the gate terminals) must follow a probability distribution typical of critical branching processes. We first note that at the initial time of the dynamics a general DMM finds itself in an unstable (unsatisfied) state. The voltages at the different terminals of the gates then start evolving, and at some intervals of time some of them undergo sudden transitions, thus creating avalanches (see typical voltage trajectories in Fig. 2, top panel) Di Ventra et al. (2017). Additionally, the memory variables in DMMs have much slower dynamics than the voltage variables Traversa and Di Ventra (2017); Di Ventra and Traversa (2018). This implies that each avalanche is independent of the others generated (mean-field condition).

Now, every time a given voltage flips from +1 to $-1$, or vice-versa, so that its corresponding Boolean variable changes from logical 1 to logical 0, or the reverse, on average, it will only have enough strength (power) to affect one other voltage in the circuit (its “offspring”). In turn, this “offspring” voltage, on average, will have enough strength to only affect at most one other voltage at the next time step, and so on.

Since all voltages in the circuit are equally important, the distribution of the number of voltages affected by a given voltage must be the same for each individual voltage at every time step (a “generation” step), and independent of both the number of flipped voltages at that time step and the number of affected voltages (offspring). Therefore, the flipping of a single voltage gives rise to a Poisson-distributed process where the average number of affected variables is $\mu\rightarrow 1$.

Under these conditions, the number of “descendants” of a flipped voltage (the size of the avalanches) is an integer random variable, $S$, described by the Borel distribution $p_{S}=(\mu S)^{S-1}e^{-\mu S}/S!$ Tanner (1961). The expectation value of $S$ is given by $\langle S\rangle=1/(1-\mu)$. Therefore, due to $\mu\rightarrow 1$, DMMs must showcase a critical branching process. In fact, in the limit of $\mu\rightarrow 1$, the Stirling approximation of the Borel distribution is proportional to $S^{-3/2}$, namely a scale-free distribution di Santo et al. (2017); Kalle Kossio et al. (2018).

III DMM equations of motion

We can now corroborate this theoretical analysis with actual numerical results. We first design a DMM for solving 3-SAT problems. Since the choice of the dynamical system representing a DMM is not unique Traversa and Di Ventra (2017); Di Ventra and Traversa (2018), we choose one very similar to the one employed in Ref. Sheldon et al. (2018) to find the ground state of the Ising spin-glasses.

In a 3-SAT problem, clauses take the form $(\ell_{i}\ \ell_{j}\ \ell_{k})$, where $\ell_{i}$ is a literal associated with a Boolean propositional variable $y_{i}$, and can be either $\ell_{i}=y_{i}$ or $\ell_{i}=\lnot y_{i}$. The variables, $y_{i}$, are transformed to continuous variables, $v_{i}$, representing terminal voltages on the self-organizing OR gates (see Fig. 1) Traversa and Di Ventra (2017); Bearden et al. (2018). The voltages are bounded, $v_{i}\in[-1,1]$, with $v_{i}\geq 0$ transformed to $y_{i}=1$ and $v_{i}<0$ transformed to $y_{i}=0$. The negation operation used by literals is trivially performed on the voltages by multiplying them by $-1$.

We then convert the $n$-th clause to a dynamical clause by interpreting literals as voltages,

[TABLE]

where subtraction is used if $\ell_{i}=y_{i}$ and addition is used if $\ell_{i}=\lnot y_{i}$, with $C_{n}\in[0,1]$. When the clause is satisfied we have $C_{n}=0$.

The dynamics of the voltages are influenced by the dynamical clauses in which the voltages appear Sheldon et al. (2018),

[TABLE]

where the sum is taken over all clauses, $C_{n}$, in which $v_{i}$ is present. The initial condition for the voltages is chosen randomly in the interval $[-1,1]$, and the solution is found when $C_{n}=0$ for all clauses. The memory variables, $x_{s,n}$ and $x_{f,n}$, along with the functions $G_{n,i}$ and $R_{n,i}$, are discussed below.

Each clause has its own memory variables, $x_{f,n}$ and $x_{s,n}$, containing information of past dynamics. The “fast” memory variable, $x_{f,n}$, determines the state of satisfiability of the clause in the recent past. Its dynamics are governed by

[TABLE]

where we have chosen $\beta=1/2$ and $\gamma=1/8$. The fast memory variable is bounded, $x_{f,n}\in[0+\epsilon,1-\epsilon]$, with the offset, $\epsilon=10^{-3}$, such that $x_{f,n}=0+\epsilon$ is interpreted to mean the clause has been satisfied for a period of time in the recent past, and $x_{f,n}=1-\epsilon$ means the clause has been unsatisfied for a period of time in the recent past. The offset is used to remove spurious steady-state solutions from Eq. (3).

The role $x_{f,n}$ has on Eq. (2) is to switch between the first and second terms in the summation. It can be seen that $x_{f,n}$ continuously switches between two modes: search for a satisfying assignment and hold the satisfying assignment. The first term in the summation contains a “gradient-like” function, $G_{n,i}(v_{i},v_{j},v_{k})$, that tries to satisfy clause $n$ by changing $v_{i}$,

[TABLE]

where the sign is chosen as in Eq. (1).

The second term in the summation of Eq. (2) contains a “rigidity” function Sheldon et al. (2018), $R_{n,i}(v_{i},v_{j},v_{k})$, which either tries to pull $v_{i}$ towards an assignment ($v_{i}=\pm 1$) that makes $C_{n}=0$ if $v_{i}$ is the voltage closest to satisfying the clause, or does nothing to $v_{i}$, if $v_{j}$ or $v_{k}$ is the voltage closest to satisfying the clause, namely

[TABLE]

Again, the signs are chosen as in Eq. (1).

The slow memory variable, $x_{s,n}$, adds weight to the gradient-like term for clause $n$ in Eq. (2), but does not affect the rigidity term. The additional weight promotes the ability to overcome the rigidity terms associated with other clauses. In essence, $x_{s,n}$ acts like a memory-assisted current generator that injects current into the circuit to push the DMM towards the solution, as originally conceived in Ref. Traversa and Di Ventra (2017). This variable is also bounded, $x_{s,n}\in[1,M]$, where $M$ is the number of variables in the problem111In practice, it is not necessary to have such a large upper bound., and we choose as its dynamics:

[TABLE]

where we have chosen $\alpha=1/100$. The slow memory variable will grow while its associated clause is unsatisfied, thereby giving the literals within that clause added weight to influence the dynamics of the voltages. In effect, $x_{s,n}$ contains memory of how often the clause was unsatisfied while traversing the phase space. We choose to initialize both the slow and fast memory variables as 1 for all clauses.

IV Results

With these DMM equations we have solved benchmark problems from previous SAT competitions. The random 3-SAT benchmarks were taken from www.satcompetition.org, and correspond to a ratio between clauses and variables of $4.3$. We have solved instances of random 3-SAT for variable sizes 420, 460, 500, 540, and 600. For each variable size we have extracted data from 100 solutions found from random initial conditions. Once a solution is found, we analyze the data for avalanches.

To identify an avalanche we analyze the magnitude of the time derivative of the vector of all the voltages, $|d\vec{v}(t)/dt|$, which we refer to as the “voltage derivative”. The random initial conditions cause a large spike in $|d\vec{v}(t)/dt|$, as seen in the bottom panel of Fig. 2. This short transient is then ignored by removing the first 5 units of time from the voltage derivative (leftmost red vertical line). When the DMM has found a solution, the voltage derivative will approach zero. We then choose to ignore also the last 10 units of time of the voltage derivative (rightmost red vertical line), so the minimum of the voltage derivative is not found at the time boundary of the simulation222Occasionally, the interval of 10 units of time is not large enough. When this occurs, the cutoff is incremented by 1 unit of time, and the minimum is checked again. The process is repeated until the minimum is not found on the time boundary..

We show in the top panel of Fig. 2 all voltages in the solution of a random 3-SAT problem. From the figure it is easy to identify “open” regions that are best described by a lack of collective flipping of voltages. These regions are characteristic of an absence of avalanche events. The identified avalanche events are enclosed in blue boxes.

Once the short initial transient and the very final approach to solution are eliminated from the voltage derivative, we find its maximum and minimum. A threshold is calculated by adding 25% of the distance between minimum and maximum to the minimum of the voltage derivative. When the voltage derivative rises above the threshold we assume that an avalanche begins, and when the voltage derivative drops below the threshold the avalanche ends. Between these two events, we check how many voltages change sign. The size of the avalanche, $S$, is defined as the number of voltages that change sign between threshold crossings.

Our scheme for determining avalanches is, of course, not without uncertainty. For instance, it is possible for multiple avalanches to be identified as a single avalanche when they are close together in time, and thresholding may also miss the voltages that flipped immediately before and after an avalanche event.

To account for the uncertainty introduced by thresholding we choose to bin the data of our numerical simulations. The bin size, $R$, is chosen using the Freedman-Diaconis rule Freedman and Diaconis (1981), $R=2IQR(S)A^{-1/3}$, rounded to the nearest integer, where $IQR(S)$ is the interquartile range of the distribution of avalanche sizes, $S$, and $A$ is the number of avalanches observed. For the 600-variable instances, the Freedman-Diaconis rule rounds to $R=9$. For better comparison, we have applied this bin size to all variable sizes.

The results of the analysis are shown in Fig. 3. Using a power-law fit we find that the initial portion of the distribution is proportional to $S^{-\tau}$, with $\tau=1.51\pm 0.02$, which is in very good agreement with the one predicted by our mean-field theory, thus giving support to the hypotheses made in that analysis.

Note that the distribution in Fig. 3 has the characteristic “bump” seen in many finite-size power-law distributions Pruessner (2012). We attribute the bump to mis-identifications made by the thresholding process. For example, in Fig. 2 we see for $100<t<125$ the thresholding has identified one avalanche. If the thresholding were to be raised slightly, the thresholding would identify three avalanches. However, we have found that changing the thresholding percentage has a negligible effect on the scale-free trend before the bump, because the composite avalanches are comprised of smaller avalanches, and thus different thresholds simply correspond to different sampling from the same distribution.

V Conclusions

In conclusion, we have provided analytical arguments, supported by numerical results, that memcomputing machines (machines that use memory to process information) undergo a critical branching process with exponent $3/2$ during their transient dynamics. The dynamics of DMMs then share some of the same features observed in many non-equilibrium phenomena encountered in Nature, and demonstrate the rich phenomenology these dynamical systems showcase, which is behind their ability to solve complex problems efficiently.

Acknowledgments – This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1650112. S.B. acknowledges support from the Alfred P. Sloan Foundation’s Minority Ph.D. Program. M.D. acknowledges partial support from the Center for Memory and Recording Research at UC San Diego.