Directed percolation and numerical stability of simulations of digital memcomputing machines
DDirected percolation and numerical stability of simulations of digital memcomputing machines
Directed percolation and numerical stability of simulations of digitalmemcomputing machines
Yuan-Hang Zhang a) and Massimiliano Di Ventra b) Department of Physics, University of California, San Diego, CA 92093, USA (Dated: 9 February 2021)
Digital memcomputing machines (DMMs) are a novel, non-Turing class of machines designed to solve combinatorialoptimization problems. They can be physically realized with continuous-time, non-quantum dynamical systems withmemory ( time non-locality ), whose ordinary differential equations (ODEs) can be numerically integrated on moderncomputers. Solutions of many hard problems have been reported by numerically integrating the ODEs of DMMs,showing substantial advantages over state-of-the-art solvers. To investigate the reasons behind the robustness and ef-fectiveness of this method, we employ three explicit integration schemes (forward Euler, trapezoid and Runge-Kutta4th order) with a constant time step, to solve 3-SAT instances with planted solutions. We show that, (i) even if most ofthe trajectories in the phase space are destroyed by numerical noise, the solution can still be achieved; (ii) the forwardEuler method, although having the largest numerical error, solves the instances in the least amount of function evalua-tions; and (iii) when increasing the integration time step, the system undergoes a “solvable-unsolvable transition” at acritical threshold, which needs to decay at most as a power law with the problem size, to control the numerical errors.To explain these results, we model the dynamical behavior of DMMs as directed percolation of the state trajectory inthe phase space in the presence of noise. This viewpoint clarifies the reasons behind their numerical robustness andprovides an analytical understanding of the unsolvable-solvable transition. These results land further support to theusefulness of DMMs in the solution of hard combinatorial optimization problems.
Certain classical dynamical systems with memory, calleddigital memcomputing machines, have been suggested tosolve a wide variety of combinatorial optimization prob-lems . In particular, simulations of their ordinary dif-ferential equations (ODEs) on traditional computers, us-ing explicit methods of integration, have shown substan-tial advantages over state-of-the-art algorithms. These re-sults are remarkable considering the fact that the ODEsof these dynamical systems are stiff, hence should becomehighly unstable if integrated with explicit methods. In thiswork, by drawing a connection between directed percola-tion and state dynamics in phase space, we provide bothnumerical and analytical reasons why explicit methods ofintegration for these ODEs are robust against the unavoid-able numerical errors. I. INTRODUCTION
MemComputing is a recently proposed (non-Turing) com-puting paradigm in which time non-locality (memory) accom-plishes both tasks of information processing and storage . Its digital (hence scalable) version (digital memcomputing ma-chines or DMMs) has been introduced specifically to solvecombinatorial optimization problems .The basic idea behind DMMs is that, instead of solvinga combinatorial optimization problem in the traditional algo-rithmic way, one maps it into a specifically designed dynam-ical system with memory, so that the only equilibrium points a) Electronic mail: [email protected] b) Electronic mail: [email protected] of that system represent the solutions of the given problem.This is accomplished by writing the problem in Boolean (oralgebraic) form, and then replacing the traditional (unidirec-tional) gates with self-organizing gates , that can satisfy theirlogical (or algebraic) relation irrespective of the direction ofthe information flow, whether from the traditional input or thetraditional output: they are terminal-agnostic gates . . Inaddition, the resulting dynamical systems can be designed sothat no periodic orbits or chaos occur during dynamics .DMMs then map a finite string of symbols into a finite string of symbols, but operate in continuous time . As a conse-quence, they seem ideally suited for a hardware implementa-tion. In fact, they can be realized in practice with non-linear, non-quantum dynamical systems with memory, such as elec-trical circuits implemented with conventional complementarymetal–oxide–semiconductor technology .On the other hand, DMMs, being classical dynamical sys-tems, are such that their state trajectory belongs to a topolog-ical manifold, known as phase space , whose dimension, D ,scales linearly with the number of degrees of freedom . Inaddition, their state dynamics are described by ordinary dif-ferential equations (ODEs) . One can then attempt to nu-merically integrate these ODEs on our traditional computers,using any integration scheme . However, naively, one wouldexpect that the computational overhead of integrating theseODEs, coupled with large numerical errors, would require anunreasonable amount of CPU time as the problem size grew,hence limiting the realization of DMMs (like quantum com-puters) to only hardware.To make things worse, the ODEs of DMMs are stiff , mean-ing that they have several, quite distinct, intrinsic time scales.This is because the memory variables of these machines havea much slower dynamics than the degrees of freedom repre-senting the logical symbols (variables) of the problem . StiffODEs are notorious for requiring implicit methods of inte- a r X i v : . [ c s . ET ] F e b irected percolation and numerical stability of simulations of digital memcomputing machines 2gration, which, in turn, require costly root-finding algorithms(such as the Newton’s method), thus making the numericalsimulations very challenging .This leads one to expect poor performance when using ex-plicit methods, such as the simplest forward Euler method,on the ODEs of the DMMs. Instead, several results using theforward Euler integration scheme on DMMs have shown thatthese machines still find solutions to a variety of combinato-rial optimization problems, including, Boolean satisfiability(SAT) , maximum satisfiability , spin glasses , machinelearning , and integer linear programming .For instance, in Ref. 7, the numerical simulations of theDMMs, using forward Euler with an adaptive time step, sub-stantially outperform traditional algorithms , in the solu-tion of 3-SAT instances with planted solutions. The resultsshow a power-law scaling in the number of integration stepsas a function of problem size, avoiding the exponential scalingseen in the performance of traditional algorithms on the sameinstance classes. Furthermore, in Ref. 7, it was found that theaverage size of the adaptive time step decayed as a power-law as the problem size grew.It is then natural to investigate how the power-law scaling ofintegration steps varies with the numerical integration schemeemployed to simulate the DMMs. In particular, by employinga constant integration time step, ∆ t , we would like to inves-tigate if the requisite time step to solve instances and controlthe numerical errors, continues to decay with a power-law, orwill require exponential decay as the problem grows.In this paper, we perform the above investigations whileattempting to solve some of the planted-solution 3-SAT in-stances used as benchmarks in Ref. 7. We will implementthree explicit integration methods (forward Euler, trapezoid,and Runge-Kutta 4th order) while numerically simulating theDMM dynamics to investigate the effect on scalability of thenumber of integration steps versus the number of variables ata given clause-to-variable ratio.Our numerical simulations indicate that, regardless of theexplicit integration method used, the ODEs of DMM are ro-bust against the numerical errors caused by the discretizationof time. The robustness of the simulations is due to the instan-tonic dynamics of DMMs , which connect critical pointsin the phase space with increasing stability. Both the criti-cal points and the instantons connecting them are of topolog-ical character . Therefore, if the integration method pre-serves critical points (in number and index) then the solutionsearch is “topologically protected” against reasonable numer-ical noise . In particular, we find that the solution can stillbe found as long as more than 1 / D of the (instantonic) pathsremains, where D is the dimension of the phase space.For each integration method, we find that when increasingthe integration time step ∆ t , the system undergoes a “solvable-unsolvable transition” at a critical ∆ t c , which decays at mostas a a power law with the problem size, avoiding the unde-sirable exponential decay that severely limits numerical sim-ulations. We also find that, even though higher-order inte-gration schemes are more favorable in most scientific prob-lems, the first-order forward Euler method works the best forthe dynamics of DMMs, providing best scalability in terms of function evaluations vs. the size of the problem. We at-tribute this to the fact that the forward Euler method preservescritical points, irrespective of the size of the integration timestep, while higher order methods may introduce “ghost criticalpoints" in the system , disrupting the instantonic dynamics ofDMMs.Finally, we provide a physical understanding of these re-sults, by showing that, in the presence of noise, the dynamicalbehavior of DMMs can be modeled as a directed percolation of the state trajectory in the phase space, with the inverse ofthe time step ∆ t playing the role of percolation probability. Wethen analytically show that, with increasing percolation prob-ability, the system undergoes a permeable-absorbing phasetransition , which resembles the “solvable-unsolvable transi-tion” in DMMs. All together, these results clarify the reasonsbehind the numerical robustness of the simulations of DMMs,and further reinforce the notion that these dynamical systemswith memory are a useful tool for the solution of hard com-binatorial optimization problems, not just in their hardwareimplementation but also in their numerical simulation.This paper is organized as follows. In Sec. II we reviewthe DMMs used in Ref. 7 to find satisfying assignments to 3-SAT instances. (Since our present paper is not a benchmarkpaper, the interested reader is directed to Ref. 7 for a com-parison of the DMM approach to traditional algorithms for3-SAT.) In Sec. III we compare the results on the scalabilityof the number of steps to reach a solution for three explicitmethods: forward Euler, trapezoid, and Runge-Kutta 4th or-der. For all three methods, our simulations show that increas-ing ∆ t (thereby increasing numerical error) induces a sort of“solvable-unsolvable phase transition,” that becomes sharperwith increasing size of the problem. However, we also showthat the “critical” time step, ∆ t c , decreases as a power law asthe size of the problem increases for a given clause-to-variableratio. In Sec. IV, we model this transition as a directed per-colation of the state trajectory in phase space, with ∆ t playingthe role of the inverse of the percolation probability. We con-clude in Sec. V with thoughts about future work. II. SOLVING 3-SAT WITH DMMS
In the remainder of this paper, we will focus on solvingsatisfiable instances of the 3-SAT problem , which is de-fined over N Boolean variables { y i = , } constrained by M clauses. Each clause consists of 3 (possibly negated) Booleanvariables connected by logical OR operations, and an instanceis solved when an assignment of { y i } is found such that all M clauses evaluate to TRUE (satisfiable).For completeness, we briefly review the dynamical systemwith memory representing our DMM. The reader is directedto Ref. 7 for a more detailed description and its mathematicalproperties.To find a satisfying assignment to a 3-SAT instance witha DMM, we transform it into a Boolean circuit, where theBoolean variables y i are transformed into continuous variablesand represented with terminal voltages v i ∈ [ − , ] (in arbi-trary units), where v i > y i =
1, and v i < y i =
0. We use ( l m , i ∨ l m , j ∨ l m , k ) to representthe m -th clause, where l m , i = y i or ¯ y i depending on whether y i is negated in this clause. Each Boolean clause is also trans-formed into a continuous constraint function: C m ( v i , v j , v k ) =
12 min (cid:0) − q m , i v i , − q m , j v j , − q m , k v k (cid:1) , (1)where q m , i = l m i = y i and q m , i = − l m i = ¯ y i . We canverify that the m -th clause evaluates to true if and only if C m < . .The idea of DMMs is to propagate information in theBoolean circuits ‘in reverse’ by using self-organizing logicgates (SOLGs) (see for an application of algebraic ones).SOLGs are designed to work bidirectionally, a property re-ferred to as “terminal agnosticism” . The terminal that is tra-ditionally an output can now receive signals like an input ter-minal, and the traditional input terminals will self-organizeto enforce the Boolean logic of the gate. This extra fea-ture requires additional (memory) degrees of freedom withineach SOLG. We introduce two additional memory variablesper gate: “short-term” memory x s , m and “long-term” mem-ory x l , m . For a 3-SAT instance, the dynamics of our self-organizing logic circuit (SOLC) are governed by Eq. (2) :˙ v i = ∑ m x l , m x s , m G m , i + (cid:0) + ζ x l , m (cid:1) ( − x s , m ) R m , i , ˙ x s , m = β ( x s , m + ε ) ( C m − γ ) , ˙ x l , m = α ( C m − δ ) , (2)where x s ∈ [ , ] , x l ∈ [ , M ] , α = , β = , γ = . , δ = . , ε = − , ζ = .
1. The “gradient-like” term ( G m , i = q m , i min ( − q m , j v j , − q m , k v k ) ) and the “rigidity” term( R m , i = ( q m , i − v i ) if C m = ( − q m , i v i ) , and R m , i =
0, other-wise) are discussed in Ref. 7 along with other details.
III. NUMERICAL SCALABILITY WITH DIFFERENTINTEGRATION SCHEMESA. Numerical details
To simulate Eq. (2) on modern computers, we need anappropriate numerical integration scheme, which necessitatesthe discretization of time. In Ref. 7 we applied the first-orderforward Euler method, which, in practice, tends to introducelarge numerical errors in the long-time integration. Yet, themassive numerical error in the trajectory of the ODE solutiondoes not translate into logical errors in the 3-SAT solution.The algorithm in Ref. 7 was capable of showing power-law scalability in the typical-case of clause distribution control(CDC) instances , though, using an adaptive time step.To further understand this robustness, we study here the be-havior of Eq. (2) under different explicit numerical integrationmethods using a constant time step, ∆ t . We will investigatethe effect of increasing the time step. Writing the ODEs as˙ x ( t ) = F ( x ( t )) , with x the collection of voltage and memory FIG. 1. Solution of clause distribution control 3-SAT instances with10 variables at clause-to-variable ratio α r = ∆ t . We performed 10 solution trials (100 instances and 100different initial conditions for each instance) for each ∆ t . We observethe number of unsolved cases decays as integration steps increase,with the solid line representing the result for all 10 trials, and shadedarea representing one standard deviation over the 100 curves calcu-lated separately using the 100 instances, with each curve covering100 different initial conditions using one specific instance. At largenumber of integration steps, the number of unsolved cases reaches aplateau, whose height only depends on the integration time step ∆ t .The fraction of solved instances represents the size of the basin ofattraction of Eqs. (2), which is similar for all instances at a certain ∆ t and shrinks as ∆ t increases. This indicates that rather than beingproblem-specific, our numerical algorithm is a general incompletesolver for 3-SAT problems. variables, and F the flow vector field , we use the followingexplicit Runge-Kutta time step : x n + = x n + ∆ t q ∑ i = ω i k i (3)where k i = F ( x n + ∆ t ∑ i − j = λ i j k j ) . Specifically, we apply threedifferent integration schemes: forward Euler method with q = ω = trapezoid method with q = ω = ( , ) , λ = method (RK4) with q = ω =( , , , ) , and λ = . (4)As an illustration, we focus on planted-solution 3-SAT in-stances at clause-to-variable ratio α r = . These instancesrequire exponential time to solve using the local-search algo-rithm walk-SAT and other state-of-the-art solvers . In theAppendix, we will show similar results for α r = N = , and solve them by numerically integratingirected percolation and numerical stability of simulations of digital memcomputing machines 4 t
0 0.51 0 0.51 0 0.51 A EulerTrapezoidRunge-Kutta 4 N=100N=200N=500N=1000N=2000N=5000N=10000
FIG. 2. The size of the basin of attraction, A , versus ∆ t , for differentnumbers of variables N and different explicit integration schemes,is well fitted by sigmoid-like curves, that become sharper as N in-creases. This indicates a “solvable-unsolvable phase transition” ata critical ∆ t c . Each data point is calculated based on 1000 solutiontrials (100 instances and 10 initial conditions per instance), and thesolid curves are fitted using the function A = / [ + exp ( − c ( ∆ t − d ))] , with c and d fitting parameters. N t c Euler t c =0.96N -0.13 Trapezoid t c =1.61N -0.14 RK4 t c =2.37N -0.13 FIG. 3. In this figure, ∆ t c is extracted from Fig. 2 where A ( ∆ t c ) = / ∆ t c shows a power-law scaling with number of variables N , and,as expected, integration methods with higher orders have a larger ∆ t c . Eqs. (2) using the forward Euler method with different val-ues of ∆ t . For each 3-SAT instance, we made 100 individualattempts with different random initial conditions, so that thetotal number of solution attempts is 10 . In Fig. 1, we see thenumber of unsolved instances decreases rapidly until reachinga plateau. When the simulations are performed again with asmaller time step, the plateau height decreases. Once reached,the plateau persists while increasing the number of integration I n t eg r a t i on s t ep s
50% solved
EulerTrapezoidRK4
90% solved
EulerTrapezoidRK4 N F un c t i on e v a l ua t i on s EulerTrapezoidRK4 N FIG. 4. The scalability curves for the three different explicit in-tegration schemes considered in this work, in which the integra-tion time step ∆ t is calculated using the power-law fit of ∆ t where A ( ∆ t ) = .
95, times an additional “safety factor” of 0.6 (see text).All three methods show a power-law scaling, and the values of thescaling exponents are 0.34 (Euler 50%), 0.50 (Trapezoid 50%), 0.44(RK4 50%), 0.24 (Euler 90%), 0.56 (Trapezoid 90%), 0.30 (RK490%). The scaling is the same, for each integration scheme, in termsof either number of steps or function evaluations; only the pre-factorchanges. For N ∈ [ , ] , Runge-Kutta 4th order (RK4) requiresthe least number of integration steps. However, in terms of functionevaluations, the simplest forward Euler method is the most efficientone. steps. B. Basin of attraction
We argue that the plateaus in Fig. 1 are caused by a re-duction of the basin of attraction in Eqs. (2), created by thenumerical integration method, rather than the hardness of theinstances themselves. To show this, first note that all solu-tion attempts succeed when ∆ t = .
15 (basin of attraction isthe entire phase space), while almost all attempts fail when ∆ t = .
35 (basin of attraction shrinks to zero size). When ∆ t falls in between, for each 3-SAT instance, the number of un-solved attempts is centered around the shaded area in Fig. 1.That is, at a certain ∆ t , the size of the basin of attraction fordifferent instances is approximately the same.We use A to denote the ratio between the volume of thebasin of attraction and the volume of the entire phase space. InFig. 2, we estimate this quantity by calculating the fraction ofsolved instances for each number of variables N and time step ∆ t for 100 different 3-SAT instances per size N , and 10 initialconditions per instance. At smaller ∆ t , all instances with alldifferent initial conditions are solved, i.e., the basin of attrac-tion is the entire phase space, A →
1. This is in agreementwith theoretical predictions for continuous-time dynamics .On the other hand, when ∆ t is too large, the structure of theirected percolation and numerical stability of simulations of digital memcomputing machines 5phase space gets modified by the numerical error introducedduring discretization, and A →
0. This is a well known effectthat has been studied also analytically in the literature .Plotting A versus ∆ t for different number of variables N anddifferent integration schemes, we get a series of sigmoid-likecurves (Fig. 2), that become sharper as N increases. Thissuggests the existence of a “phase transition”: when goingbeyond a critical ∆ t c , there is a drastic change in A , and thesystem undergoes a solvable-unsolvable transition . C. Integration step scalability
The above results allow us to determine how the integrationstep, ∆ t , for the three integration schemes we have chosen,scales as the size of the problem increases. In Fig. 3, wedefine ∆ t c such that A ( ∆ t c ) = /
2, and determine the relationbetween ∆ t c and N .We find ∆ t c scales as a power law with the number of vari-ables N . This is a major result because it shows that we do notneed ∆ t to decrease exponentially to have the same successrate. In Ref. 7, it was demonstrated, using topological fieldtheory , that in the ideal (noise-free) continuous-time dy-namics, the physical time required for our dynamical systemto reach the solution scales as O ( N γ ) , with γ ≤ ∆ t c ∼ N − δ . Coupled with the previous analytical re-sults , this means that the number of integration steps scalesas O ( N γ + δ ) . In other words, discretization only adds a O ( N δ ) overhead to the complexity of our algorithm, indicating thatDMMs can be efficiently simulated on classical computers .In Fig. 3, note that ∆ t c > − for all three integration meth-ods. This is quite unexpected for a stiff system simulated withexplicit integration methods, because a large time step intro-duces large local truncation errors at each step of the integra-tion. This error accumulates and should destroy the trajectory of the ODEs we are trying to simulate. However, we can stillsolve all tested instances with Eqs. (2) even at such a large ∆ t . In Sec. IV, we will provide an explanation of why this ispossible.Here, we note that choosing an appropriate ∆ t could speedup the algorithm significantly, as smaller ∆ t leads to excessiveintegration steps, while larger ∆ t may cause the solver to fail.To find an appropriate time step for our scalability tests, wechoose ∆ t from Fig. 2 such that 95% of the trials have founda solution, and multiply that ∆ t by a “safety factor” of 0 . A ( ∆ t ) ≈ ∆ t , we plot in Fig. 4the typical-case and 90th percentile scalability curves up to N = . For each variable size, N , we perform 1000 3-SATsolution trials (100 instances and 10 initial conditions each),and report the number of integration steps when 50% and 90%of the trials are solved. Additionally, we report the number of function evaluations (namely, the number of times the right-hand side of Eqs. (2) is evaluated), which differs from thenumber of integration steps by only a constant.Both quantities show a power-law scalability for all three integration schemes. Since ∆ t c is larger for higher-orderedODE solvers, the latter ones typically require fewer steps tosolve the problems. However, when taking the total number offunction evaluations into account, the simplest forward Eulermethod becomes the most efficient integration scheme for ourODEs. IV. DIRECTED PERCOLATION AND NOISE IN DMMS
As anticipated, our results are unexpected in view of thelarge truncation errors introduced during the simulations ofthe ODEs of the DMMs. However, one might have predictedthe numerical robustness upon considering the type of dynam-ics the DMM executes in phase space.
A. Instantonic dynamics
It was shown in Refs. , using (supersymmetric) topo-logical field theory, that the only “low-energy”, “long-wavelength” dynamics of DMMs is a collection of elementary instantons (a composite instanton). The latter ones are topo-logically non-trivial classical trajectories that connect criticalpoints (those for which the flow vector field is equal to zero)with decreasing index (number of unstable directions) .Critical points are also of topological character. For instance,their number cannot change unless we change the topology ofphase space .In addition, given two distinct (in terms of indexes) criti-cal points, there are several (a family of) trajectories (instan-tons) that may connect them, since the unstable manifold ofthe “initial” critical point may intersect the stable manifold ofthe “final” critical point at several points in the phase space.In the ideal (noise-free) continuous time dynamics, if the onlyequilibria (fixed points of the dynamics) are the solutions ofthe given problem (as shown, e.g., in ), then the state trajec-tory is driven towards the equilibria by the voltages that setthe input of the problem the DMM attempts to solve. B. Physical noise
In the presence of (physical) noise instead, anti-instantons are generated in the system. These are (time-reversed) trajec-tories that connect critical points of increasing index. How-ever, anti-instantons are gapped , namely their amplitude is ex-ponentially suppressed with respect to the corresponding am-plitude of the instantons . This means that the “low-energy”,“long wave-length” dynamics of DMMs is still a successionof instantons, even in the presence of (moderate) noise.Nevertheless, suppose an instanton connects two criticalpoints, and immediately after an anti-instanton occurs for the same critical points (even if the two trajectories are different).For all practical purposes, the “final” critical point of the orig-inal instanton has never been reached, and that critical pointcan be called absorbing , in the sense that the trajectory would“get stuck” on that one, or wanders in other regions of theirected percolation and numerical stability of simulations of digital memcomputing machines 6phase space. In other words, in the presence of noise, there isa finite probability for some state trajectory to explore a muchlarger region of the phase space. The system then needs toexplore some other (anti-)instantonic trajectory to reach theequilibria, solutions of the given problem. Nevertheless, dueto the fact that anti-instantons are gapped, and the topologi-cal character of both critical points and instantons connectingthem, if the physical noise level is not high enough to changethe topology of the phase space, the dynamical system wouldstill reach the equilibria. C. Directed percolation
If we visualize the state trajectory as the one traced by aliquid in a corrugated landscape, the above suggests an in-triguing analogy with the phenomenon of directed percolation (DP) , with the critical points acting as ‘pores’, and instan-tons as ‘channels’ connecting ‘neighboring’ (in the sense ofindex difference) critical points.DP is a well-studied model of a non-equilibrium (continu-ous) phase transition from a fluctuating permeable phase to anabsorbing phase . It can be intuitively understood as a liq-uid passing through a porous substance under the influence ofa field (e.g., gravity). The field restricts the direction of theliquid’s movement, hence the term “directed”. In this model,neighboring pores are connected by channels with probability p , and disconnected otherwise. When increasing p , the sys-tem goes through a phase transition from the absorbing phaseinto a permeable phase at a critical threshold p c . D. Numerical noise
In numerical ODE solvers, discrete time also introducessome type of noise (truncation errors), and the considerationswe have made above on the presence of anti-instantons stillhold. However, numerical noise could be substantially moredamaging than physical noise. The reason is twofold.First, as we have already anticipated, numerical noise accu-mulates during the integration of the ODEs. In some sense, itis then always non-local . Physical noise, instead, is typicallylocal (in space and time), hence it is a local perturbation of thedynamics . As such, if it is not too large, it cannot change thephase space topology.Second, unlike physical noise, integration schemes maychange the topology of phase space explicitly . This is because,when one transforms the continuous-time ODEs (2) into theirdiscrete version (a map ), this transformation can introduce ad-ditional critical points in the phase space of the map, whichwere not present in the original phase space of the continu-ous dynamics. These extra (undesirable) critical points aresometimes called ghosts . For instance, while the forwardEuler method can never introduce such ghost critical points,irrespective of the size of ∆ t (because both the continuous dy-namics and the associated map have the same flow field F ),both the trapezoid and Runge-Kutta 4th order may do so if ∆ t is large enough. These critical points would then furtherdegrade the dynamics of the system.With these preliminaries in mind, let us now try to quantifythe analogy between the DMM dynamics in the presence ofnumerical noise and directed percolation. E. Paths to solution
First of all, the integration step, ∆ t , must be inversely re-lated to the percolation probability p : when ∆ t tends to zero,the discrete dynamics approach the ideal (noise-free) dynam-ics for which p →
1. In the opposite case, when ∆ t increases, p must decrease.In directed percolation, the order parameter is the densityof active sites . This would correspond to the density of“achievable” critical points in DMMs. However, due to thevast dimensionality of their phase space (even for relativelysmall problem instances), this quantity is not directly accessi-ble in DMMs.Instead, what we can easily evaluate is the ratio of suc-cessful solution attempts: starting from a random point in thephase space (initial condition of the ODEs (2)), and letting thesystem evolve for a sufficiently long time, a path would eithersuccessfully reach the solution (corresponding to a permeablepath in DP) or fail to converge (corresponding to an absorb-ing path in DP). By repeating this numerical experiment fora large number of trials, the starting points of the permeablepaths essentially fill the basin of attraction of our dynamicalsystem. The advantage of considering this quantity is that theratio of permeable paths in bond DP models can be calculatedanalytically.Consider then a 3-SAT problem with only one solution. ADMM for such a problem can reach the solution from severalinitial conditions. This translates into a D -dimensional cone-shaped square lattice for the DP model (see Fig. 5). The (top)base of the cone would represent the starting points, and theapex of the cone would represent the solution point. A per-meable path connects the base to the apex, and an absorbingpath ends in the middle of the lattice with all bonds beneath itdisconnected.To begin with, we assume that all starting points are occu-pied with equal probability, and calculate the expectation ofthe number of permeable and absorbing paths. This can bedone iteratively: the expected number of permeable paths ata given site equals the sum of permeable paths at its prede-cessor sites, times p , the percolation probability. Assumingthe number of time steps required to reach the apex is T , thenthe expected number of permeable paths (cid:104) n p , T (cid:105) = ( Dp ) T . Anillustration of the calculation on a cone-shaped 2d lattice isshown in Fig. 5.Instead, the number of absorbing paths at each site is thenumber of permeable paths at that site, times ( − p ) D , theprobability that all bonds beneath it are disconnected. (Notethat this probability has a different expression at the boundaryof the lattice, but as D and T are usually very large in DMMs,we ignore the boundary effect here.) Then, the total numberirected percolation and numerical stability of simulations of digital memcomputing machines 7 D i r ec ti on o f p e r c o l a ti on Number of permeable paths from the (top) base to the current site:
Actual number (Expectation) (1) (1) (1) (1) (1) (2p) (2p) (2p) (2p) (4p ) (4p ) (4p ) (8p ) (8p ) (16p ) FIG. 5. Illustration of permeable and absorbing paths in DP on acone-shaped 2d lattice. Connected bonds are represented with solidlines, and disconnected ones are represented with dotted lines. Apermeable path connects the base to the apex, and an absorbing pathends in the middle of the lattice with all bonds beneath it discon-nected. The number of permeable paths n p is calculated iteratively: n p at a given site equals the sum of n p at its connected predecessorsites, and the expectation (cid:104) n p (cid:105) at a given site equals the sum of (cid:104) n p (cid:105) at all its predecessor sites times p , the percolation probability. FIG. 6. The ratio of permeable paths, calculated with Eq. (9), with δ = Dp −
1, plotted for different values of the dimension D . The insetshows the same curves in log scale. The curves exhibit a sigmoidalbehavior, which is similar to the one in Fig. 2. However, note that,since ∆ t in DMMs is inversely related to the probability p in DP, theyare curving in opposite directions. of absorbing paths is (cid:104) n a , T (cid:105) = ( − p ) D T − ∑ i = ( T + − i )( Dp ) i = ( − p ) D ( − Dp ) (cid:0) ( Dp ) T ( Dp − ) + + T − Dp ( + T ) (cid:1) (5)The ratio of permeable paths r is simply (cid:104) n p , T (cid:105) / ( (cid:104) n p , T (cid:105) + (cid:104) n a , T (cid:105) ) : r = + ( − p ) D ( − Dp ) ( Dp − + ( Dp ) − T ( + T − Dp ( + T ))) . (6)Note that the transition occurs near p c = D . Therefore, letus define p ≡ + δ D , with δ ∈ R , and use δ as the new variable.To further simplify Eq. (6), we look at the limits D → ∞ and T → ∞ . Under this approximation, we have ( − p ) D = ( − + δ D ) D ≈ e − − δ , (7)and ( Dp ) − T = ( + δ ) − T ≈ (cid:26) , δ (cid:29) T ∞ , δ (cid:28) − T (8)Using Eqs. (7) and (8), Eq. (6) becomes r = (cid:40) + e − − δ ( + δ ) / δ , δ (cid:29) T , δ (cid:28) − T (9)Near δ =
0, we havelim δ → + r ( δ ) = lim δ → − r ( δ ) = δ = Dp − D . Comparing toFig. 2, we can already see a few similarities: both figuresexhibit a sigmoidal behavior. However, as ∆ t in DMMs isinversely related to p in DP, they are curving in opposite di-rections.Recall that the size of the basin of attraction, A , in DMMscorresponds to the ratio r in DP. To model their relation, wethen use the ansatz p = N (cid:0) a ∆ t − b (cid:1) , with a and b some realnumbers, and fit A to ∆ t using Eq. (9). (Note that this trialfunction has a meaning only near ∆ t c .)Fig. 7 shows the fitting result, where the number of vari-ables N ranges between 10 and 10 , and each data pointis obtained by numerically integrating Eqs. (2) for 1000 3-SAT solution attempts (100 instances and 10 initial conditionseach), until a plateau, as in Fig. 1, has been reached. Notethat the horizontal axis represents 1 / ( N ∆ t ) . The fitting re-sult for forward Euler and trapezoid are similar (see Fig. 8).The effective percolation probability p = N (cid:0) a ∆ t − b (cid:1) , where a ∼ , b ∼ N ζ , and the fitted effective dimensionality D ∼ N .The results are a bit different for RK4: The fitted dimen-sionality D becomes sub-linear in N , with D = . N . .Meanwhile, a starts to scale up with N , and the best fit wefound is a power law fit between log a and log N , with log a = . × − ( log ( N )) . . This discrepancy between RK4 andthe forward Euler/trapezoid schemes can be attributed to theproliferation of ghost critical points in the phase space ofthe RK4 map. Similarly, the reason why the forward Eulerirected percolation and numerical stability of simulations of digital memcomputing machines 8 -4 -3 -2 -1
0 0.51 0 0.51 0 0.51 A EulerTrapezoidRunge-Kutta 4
N=100N=200N=500N=1000N=2000N=5000N=10000
FIG. 7. Relation, close to ∆ t c , between the size of the basin of attrac-tion A and the discretization time step ∆ t , plotted for different systemsizes N , for the three explicit methods used in this work. Each datapoint is obtained by numerically integrating Eqs. (2) over 1000 3SATsolution attempts (100 instances and 10 initial conditions each), untilreaching a plateau. The solid curves are fitted using Eq. (6), with A corresponding to r and p = N (cid:0) a ∆ t − b (cid:1) , with a and b fitting parame-ters. method has the smallest scaling exponent of the three meth-ods (as shown in Fig. 4) may be attributed to the fact thatforward Euler cannot introduce ghost critical points.The dimensionality D is proportional to the dimensionalityof the composite instanton (namely the number of elementaryinstantons to reach the solution). In Ref. 19, this dimension-ality was predicted to scale (sub-)linearly with the number ofvariables N (in the noise-free case), and it is smaller than thedimensionality of the actual phase space ( ( + × α r ) N = N in the present case). This prediction agrees with our numericalresults.Finally, using the correspondence between DMMs in thepresence of noise and DP, we can qualitatively explain whythe DMMs numerically integrated still provide solution tothe problem they are designed to solve even with such largenumerical errors introduced during integration. In DMMs,the dimensionality of the composite instanton, D , is usuallyvery large, and DP tells us that the percolation threshold p c = / D ∼ / N .Therefore, even if most of the paths in the phase space aredestroyed by noise, we can still achieve the solution as longas more than 1 / D ∼ / N of the (instantonic) paths remain.This argument, in tandem with the fact that critical points andinstantons have a topological character , ensures the ro-bustness of DMMs in numerical simulations. V. CONCLUSIONS
In conclusion, we have shown that the dynamics of digitalmemcomputing machines (DMMs) under discrete numerical
FIG. 8. The fitted parameters a , b , and D for different variable size N for the trial function p = N (cid:0) a ∆ t − b (cid:1) , which connects the time step ∆ t to the percolation probability p . D is the dimensionality of theDP lattice, which corresponds to the dimensionality of the compositeinstanton. solvers can be described as a directed percolation of state tra-jectory in the phase space. The inverse time step, 1 / ∆ t , playsthe role of percolation probability, and the system undergoes a“solvable-unsolvable phase transition” at a critical ∆ t c , whichscales as a power law with problem size. In other words, forthe problem instances considered, we have numerically foundthat the integration time step does not need to decrease expo-nentially with the size of the problem, in order to control thenumerical errors.This result is quite remarkable considering the fact that wehave only employed explicit methods of integration (forwardEuler, trapezoid, and Runge-Kutta 4th order) for stiff ODEs.(In fact, the forward Euler method, although having the largestnumerical error, solves the instances in the least amount offunction evaluations.) It can be ultimately traced to the typeof dynamics the DMMs perform during the solution search,which is a composite instanton in phase space. Since instan-tons, and the critical points they connect, are of topologicalcharacter, perturbations to the actual trajectory in phase spacedo not have the same detrimental effect as the changes oftopology of the phase space.Since numerical noise is typically far worse than physicalnoise, these results further reinforce the notion that these ma-chines, if built in hardware, would be topologically protectedagainst moderate physical noise and perturbations.However, let us note that we did not prove that the dynamicsof DMMs with numerical (or physical) noise belong to the DPuniversality class. In fact, according to the DP-conjecture ,a given model belongs to such a universality class if (i) themodel displays a continuous phase transition from a fluctuat-irected percolation and numerical stability of simulations of digital memcomputing machines 9ing active phase into a unique absorbing state; (ii) the transi-tion is characterized by a non-negative one-component orderparameter; (iii) the dynamic rules are short-ranged; (iv) thesystem has no special attributes such as unconventional sym-metries, conservation laws, or quenched randomness.It is easy to verify that DMMs under numerical simula-tions satisfy properties (i) , (iii) , and (iv) . However, verifyingproperty (ii) is not trivial, in view of the vast phase space ofDMMs. Still, the similarities we have outlined in this paper,between DMMs in the presence of noise and DP, help us betterunderstand how these dynamical systems with memory work,and why their simulations are robust against the unavoidablenumerical errors. DATA AVAILABILITY
The 3-SAT instances and raw data used to generate all fig-ures in this paper are available upon request from the authors.
ACKNOWLEDGMENTS
We thank Sean Bearden for providing the MATLAB codeto solve the 3-SAT instances with DMMs, as well as insight-ful discussions. Work supported by NSF under grant No.2034558. All memcomputing simulations reported in this pa-per have been done on a single core of an AMD EPYC server.irected percolation and numerical stability of simulations of digital memcomputing machines 10
Appendix: Results for α r = Here, we show that the results presented in the main texthold also for other clause-to-variable ratios. As examples, wechoose α r =
6. (Similar scalability results have been reportedin Ref. 7 for α r = . ∆ t c for all clause-to-variable ratios, with ofcourse, different power laws.Figure A1 have been obtained as discussed in the main text,and show ∆ t c vs. number of variables for α r = N t c r =6 Euler t c =1.17N -0.21 Trapezoid t c =3.24N -0.30 RK4 t c =3.40N -0.24 FIG. A1. Critical ∆ t c defined as A ( ∆ t c ) = / α r = ∆ t c shows a power-law scal-ing with N , and, as expected, integration methods with higher ordershave a larger ∆ t c . In Fig. A2, instead we show the scalability curves for α r =
6, considering all three explicit integration methods. As in themain text, we find that the forward Euler method, althoughhaving the largest numerical error, solves the instances in theleast amount of function evaluations for α r =
6. Every datapoint in the curves is obtained with 100 3-SAT instances, with10 solution trials for each instance. M. Di Ventra and F. Traversa, J. Appl. Phys. , 180901 (2018). M. Di Ventra and Y. V. Pershin, Nature Physics , 200 (2013). F. L. Traversa and M. Di Ventra, Chaos: An Interdisciplinary Journal ofNonlinear Science , 023107 (2017). F. L. Traversa and M. Di Ventra, arXiv:1808.09999 (2018). S. Bearden, F. Sheldon, and M. Di Ventra, Europhysics Letters , 30005(2019). M. Di Ventra and F. L. Traversa, Phys. Lett. A , 3255 (2017). S. R. Bearden, Y. R. Pei, and M. Di Ventra, Scientific Reports , 1 (2020). O. Goldreich,
Foundations of Cryptography (Cambridge University Press,2001). T. Sauer,
Numerical Analysis (Pearson, 2018). F. L. Traversa, P. Cicotti, F. Sheldon, and M. Di Ventra, Complexity ,7982851 (2018). F. Sheldon, P. Cicotti, F. L. Traversa, and M. Di Ventra, IEEE Trans. NeuralNetw. Learn. Syst. , 2222 (2020). F. Sheldon, F. L. Traversa, and M. Di Ventra, Physical Review E ,053311 (2019). H. Manukian, F. L. Traversa, and M. Di Ventra, Neural Networks , 1(2019). I n t eg r a t i on s t ep s
50% solved
EulerTrapezoidRK4
90% solved
EulerTrapezoidRK4 N F un c t i on e v a l ua t i on s EulerTrapezoidRK4 N r =6 FIG. A2. The scalability curves at α r = ∆ t is calculated using the power-law fitof ∆ t where A ( ∆ t ) = .
95, times an additional “safety factor” of 0.6(see main text). All three methods show a power-law scaling. Thevalues of the scaling exponents are 0.62 (Euler 50%), 0.84 (Trape-zoid 50%), 0.73 (RK4 50%), 0.54 (Euler 90%), 0.77 (Trapezoid90%), 0.74 (RK4 90%). The scaling is the same, for each integrationscheme, in terms of either number of steps or function evaluations;only the pre-factor changes. For N ∈ [ , ] , Runge-Kutta 4th or-der (RK4) requires the least number of integration steps. However,in terms of function evaluations, the simplest forward Euler methodis the most efficient one. H. Manukian, Y. R. Pei, S. R. B. Bearden, and M. Di Ventra, Comm. Phys. , 1 (2020). B. Selman and H. Kautz, in
IJCAI , Vol. 93 (Citeseer, 1993) pp. 290–295. M. Mézard, G. Parisi, and R. Zecchina, Science , 812 (2002). A. Biere, Proceedings of SAT Competition 2017: Solver and BenchmarksDescriptions , 14 (2017). M. Di Ventra, F. L. Traversa, and I. Ovchinnikov, Ann. Phys. (Berlin) ,1700123 (2017). M. Di Ventra and I. V. Ovchinnikov, Annals of Physics , 167935 (2019). R. Rajaraman,
Solitons and Instantons: An Introduction to Solitons andInstantons in Quantum Field Theory (North Holland, 1982). S. Coleman,
Aspects of Symmetry, Chapter 7 (Cambridge University Press,1977). Incidentally, this is also the reason why the hardware implementation ofDMMs would be robust against reasonable physical noise . Note that all instances can be solved, as they have planted solutions. J. H. Cartwright and O. Piro, International Journal of Bifurcation and Chaos , 427 (1992). W. Barthel, A. K. Hartmann, M. Leone, F. Ricci-Tersenghi, M. Weigt, andR. Zecchina, Physical review letters , 188701 (2002). Using the method of with p = . A. M. Stuart and A. Humphries, Acta numerica , 467 (1994). Of course, this is not an analytical proof of the efficiency of the numerical simulations, just an empirical result. A. T. Fomenko,
Visual Geometry and Topology (Springer-Verlag, 1994). I. V. Ovchinnikov, Chaos , 013108 (2013). In fact, moderate noise may even help accelerate the time to solution, byreducing the time spent on the critical points’s stable directions. M. Henkel, H. Hinrichsen, S. Lübeck, and M. Pleimling,
Non-equilibriumphase transitions , Vol. 1 (Springer, 2008). N. van Kampen,
Stochastic Processes in Physics and Chemistry (Elsevier, irected percolation and numerical stability of simulations of digital memcomputing machines 11 H.-K. Janssen, Zeitschrift für Physik B Condensed Matter , 151 (1981). P. Grassberger, Zeitschrift für Physik B Condensed Matter47