Directed percolation and numerical stability of simulations of digital memcomputing machines

MemComputing is a recently proposed (non-Turing) computing paradigm in which time non-locality (memory) accomplishes both tasks of information processing and storage.² Its digital (hence scalable) version (digital memcomputing machines or DMMs) has been introduced specifically to solve combinatorial optimization problems.^1,3

The basic idea behind DMMs is that, instead of solving a combinatorial optimization problem in the traditional algorithmic way, one maps it into a specifically designed dynamical system with memory so that the only equilibrium points of that system represent the solutions of the given problem. This is accomplished by writing the problem in Boolean (or algebraic) form and then replacing the traditional (unidirectional) gates with self-organizing gates that can satisfy their logical (or algebraic) relation irrespective of the direction of the information flow, whether from the traditional input or the traditional output: they are terminal-agnostic gates.^1,3–5 In addition, the resulting dynamical systems can be designed so that no periodic orbits or chaos occur during dynamics.^6,7

DMMs then map a finite string of symbols into a finite string of symbols, but operate in continuous time, hence they are distinct from Turing machines that operate in discrete time. (Note that universal memcomputing machines have been shown to be Turing-complete, but not Turing-equivalent.⁸) As a consequence, they seem ideally suited for a hardware implementation. In fact, they can be realized in practice with non-linear, non-quantum dynamical systems with memory, such as electrical circuits implemented with conventional complementary metal–oxide–semiconductor technology.³ 

On the other hand, DMMs, being classical dynamical systems, are such that their state trajectory belongs to a topological manifold, known as phase space, whose dimension, $D$⁠, scales linearly with the number of degrees of freedom.⁹ In addition, their state dynamics are described by ordinary differential equations (ODEs).^1,3 One can then attempt to numerically integrate these ODEs on our traditional computers, using any integration scheme.¹⁰ However, naively, one would expect that the computational overhead of integrating these ODEs, coupled with large numerical errors, would require an unreasonable amount of CPU time as the problem size grew, hence limiting the realization of DMMs (like quantum computers) to only hardware.

To make things worse, the ODEs of DMMs are stiff,¹⁰ meaning that they have several, quite distinct, intrinsic time scales. This is because the memory variables of these machines have a much slower dynamics than the degrees of freedom representing the logical symbols (variables) of the problem.^1,3 Stiff ODEs are notorious for requiring implicit methods of integration, which, in turn, require costly root-finding algorithms (such as Newton’s method), thus making the numerical simulations very challenging.¹⁰ 

This leads one to expect poor performance when using explicit methods, such as the simplest forward Euler method, on the ODEs of the DMMs. Instead, several results using the forward Euler integration scheme on DMMs have shown that these machines still find solutions to a variety of combinatorial optimization problems, including Boolean satisfiability (SAT),⁷ maximum satisfiability,^11,12 spin glasses,¹³ machine learning,^14,15 and integer linear programming.⁴ 

For instance, in Ref. 7, the numerical simulations of the DMMs, using forward Euler with an adaptive time step, substantially outperform traditional algorithms,^16–18 in the solution of 3-SAT instances with planted solutions. The results show a power-law scaling in the number of integration steps as a function of problem size, avoiding the exponential scaling seen in the performance of traditional algorithms on the same instance classes. Furthermore, in Ref. 7, it was found that the average 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 of integration steps varies with the numerical integration scheme employed to simulate the DMMs. In particular, by employing a constant integration time step, $Δt$⁠, we would like to investigate if the requisite time step to solve instances and control the numerical errors continues to decay with a power law or will require exponential decay as the problem grows.

In this paper, we perform the above investigations while attempting to solve some of the planted-solution 3-SAT instances used as benchmarks in Ref. 7. We will implement three explicit integration methods (forward Euler, trapezoid, and Runge–Kutta fourth order) while numerically simulating the DMM dynamics to investigate the effect on scalability of the number of integration steps vs the number of variables at a given clause-to-variable ratio.

Our numerical simulations indicate that, regardless of the explicit integration method used, the ODEs of DMM are robust against the numerical errors caused by the discretization of time. The robustness of the simulations is due to the instantonic dynamics of DMMs,^19,20 which connect critical points in the phase space with increasing stability. Instantons in dissipative systems are specific heteroclinic orbits connecting two distinct (in index) critical points. Both the critical points and the instantons connecting them are of topological character.^21,22 Therefore, if the integration method preserves critical points (in number and index), then the solution search is “topologically protected” against reasonable numerical noise.²³ 

For each integration method, we find that when increasing the integration time step $Δt$⁠, the system undergoes a “solvable–unsolvable transition”²⁴ at a critical $Δtc$⁠, which decays at most as a a power law with the problem size, avoiding the undesirable exponential decay that severely limits numerical simulations. We also find that, even though higher-order integration schemes are more favorable in most scientific problems, the first-order forward Euler method works the best for the dynamics of DMMs, providing best scalability in terms of function evaluations vs the size of the problem. We attribute this to the fact that the forward Euler method preserves critical points, irrespective of the size of the integration time step, while higher order methods may introduce “ghost critical points” in the system,²⁵ disrupting the instantonic dynamics of DMMs.

Finally, we provide a physical understanding of these results, by showing that, in the presence of noise, the dynamical behavior of DMMs can be modeled as a directed percolation of the state trajectory in the phase space, with the inverse of the time step $Δt$ playing the role of percolation probability. We then analytically show that, with increasing percolation probability, the system undergoes a permeable-absorbing phase transition, which resembles the “solvable–unsolvable transition” in DMMs. All together, these results clarify the reasons behind the numerical robustness of the simulations of DMMs and further reinforce the notion that these dynamical systems with memory are a useful tool for the solution of hard combinatorial optimization problems, not just in their hardware implementation but also in their numerical simulation.

This paper is organized as follows. In Sec. II, we review the DMMs used in Ref. 7 to find satisfying assignments to 3-SAT instances. (Since our present paper is not a benchmark paper, the interested reader is directed to Ref. 7 for a comparison of the DMM approach to traditional algorithms for 3-SAT.) In Sec. III, we compare the results on the scalability of the number of steps to reach a solution for three explicit methods: forward Euler, trapezoid, and Runge–Kutta fourth order. For all three methods, our simulations show that increasing $Δt$ (thereby increasing numerical error) induces a sort of “solvable–unsolvable phase transition” that becomes sharper with increasing size of the problem. However, we also show that the “critical” time step, $Δtc$⁠, decreases as a power law as the size of the problem increases for a given clause-to-variable ratio. In Sec. IV, we model this transition as a directed percolation of the state trajectory in phase space, with $Δt$ playing the role of the inverse of the percolation probability. We conclude in Sec. V with thoughts about future work.

In the remainder of this paper, we will focus on solving satisfiable instances of the 3-SAT problem,²⁶ which is defined over $N$ Boolean variables ${yi=0,1}$ constrained by $M$ clauses. Each clause consists of three (possibly negated) Boolean variables connected by logical OR operations, and an instance is solved when an assignment of ${yi}$ is found such that all $M$ clauses evaluate to TRUE (satisfiable).

For completeness, we briefly review the dynamical system with memory representing our DMM. The reader is directed to Ref. 7 for a more detailed description and its mathematical properties.

To find a satisfying assignment to a 3-SAT instance with a DMM, we transform it into a Boolean circuit, where the Boolean variables $yi$ are transformed into continuous variables and represented with terminal voltages $vi∈[−1,1]$ (in arbitrary units), where $vi>0$ corresponds to $yi=1$ and $vi<0$ corresponds to $yi=0$⁠. We use $(lm,i∨lm,j∨lm,k)$ to represent the $m$th clause, where $lm,i=yi$ or $y¯i$ depending on whether $yi$ is negated in this clause. Each Boolean clause is also transformed into a continuous constraint function,

where $qm,i=1$ if $lmi=yi$ and $qm,i=−1$ if $lmi=y¯i$⁠. We can verify that the $m$th clause evaluates to true if and only if $Cm<0.5$⁠.⁷ 

The idea of DMMs is to propagate information in the Boolean circuits “in reverse” by using self-organizing logic gates (SOLGs)^1,3 (see Ref. 4 for an application of algebraic ones).

SOLGs are designed to work bidirectionally, a property referred to as “terminal agnosticism.”⁵ The terminal that is traditionally an output can now receive signals like an input terminal, and the traditional input terminals will self-organize to enforce the Boolean logic of the gate. This extra feature requires additional (memory) degrees of freedom within each SOLG. We introduce two additional memory variables per gate: “short-term” memory $xs,m$ and “long-term” memory $xl,m$⁠. For a 3-SAT instance, the dynamics of our self-organizing logic circuit (SOLC) are governed by⁷ 

where $xs∈[0,1]$⁠, $xl∈[1,104M]$⁠, $α=5,β=20,γ=0.25,δ=0.05,ϵ=10−3,andζ=0.1$⁠. The “gradient-like” term [$Gm,i=12qm,imin(1−qm,jvj,1−qm,kvk)$] and the “rigidity” term [$Rm,i=12(qm,i−vi)$ if $Cm=12(1−qm,ivi)$⁠, and $Rm,i=0$⁠, otherwise] are discussed in Ref. 7 along with other details.

To simulate Eq. (2) on modern computers, we need an appropriate numerical integration scheme, which necessitates the discretization of time. In Ref. 7, we applied the first-order forward Euler method, which, in practice, tends to introduce large numerical errors in the long-time integration. Yet, the massive numerical error in the trajectory of the ODE solution does 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 behavior of Eq. (2) under different explicit numerical integration methods using a constant time step, $Δt$⁠. We will investigate the effect of increasing the time step. Writing the ODEs as $x˙(t)=F(x(t))$⁠, with $x$ being the collection of voltage and memory variables and $F$ being the flow vector field,¹ we use the following explicit Runge–Kutta time step:¹⁰ 

where $ki=F(xn+Δt∑j=1i−1λijkj)$⁠. Specifically, we apply three different integration schemes: forward Euler method with $q=1$⁠, $ω1=1$⁠; trapezoid method with $q=2$⁠, $ω=(12,12)$⁠, $λ21=1$⁠; fourth order Runge–Kutta method (RK4) with $q=4$⁠, $ω=(16,13,13,16)$⁠, and

Note that Eq. (2) is stiff, meaning that the explicit integration methods we consider here should diverge quite fast, within a few integration steps, irrespective of the problem instances we choose.¹⁰ However, as we will show below, the explicit integration methods actually work very well for these equations, and the integration time step can even be chosen quite large.

As an illustration, we focus on planted-solution 3-SAT instances at clause-to-variable ratio $αr=8$⁠.²⁷ These instances require exponential time to solve using the local-search algorithm walk-SAT²⁶ and other state-of-the-art solvers.⁷ In  Appendix A, we will show similar results for $αr=6$⁠.

We start with 100 CDC 3-SAT instances with the number of variables $N=104$ and solve them by numerically integrating Eq. (2) using the forward Euler method with different values of $Δt$⁠. For each 3-SAT instance, we made 100 individual attempts with different random initial conditions so that the total number of solution attempts is $104$⁠. In Fig. 1, we see the number of unsolved instances decreases rapidly until reaching a plateau. When the simulations are performed again with a smaller time step, the plateau height decreases. Once reached, the plateau persists while increasing the number of integration steps.

We argue that the plateaus in Fig. 1 are caused by a reduction in the basin of attraction in Eq. (2), created by the numerical integration method, rather than the hardness of the instances themselves. To show this, first note that all solution attempts succeed when $Δt=0.15$ (basin of attraction is the entire phase space), while almost all attempts fail when $Δt=0.35$ (basin of attraction shrinks to zero size). When $Δt$ falls in between, for each 3-SAT instance, the number of unsolved attempts is centered around the shaded area in Fig. 1. That is, at a certain $Δt$⁠, the size of the basin of attraction for different instances is approximately the same.

We use $A$ to denote the ratio between the volume of the basin of attraction and the volume of the entire phase space. In Fig. 2, we estimate this quantity by calculating the fraction of solved instances for each number of variables $N$ and time step $Δt$ for 100 different 3-SAT instances per size $N$ and 10 initial conditions per instance. At smaller $Δt$⁠, all instances with all different initial conditions are solved, i.e., the basin of attraction is the entire phase space, $A→1$⁠. This is in agreement with theoretical predictions for continuous-time dynamics.⁷ On the other hand, when $Δt$ is too large, the structure of the phase space gets modified by the numerical error introduced during discretization, and $A→0$⁠. This is a well-known effect that has been studied also analytically in the literature.²⁸ 

Plotting $A$ vs $Δt$ for different number of variables $N$ and different integration schemes, we get a series of sigmoid-like curves (Fig. 2) that become sharper as $N$ increases. This suggests the existence of a “phase transition”: when going beyond a critical $Δtc$⁠, there is a drastic change in $A$⁠, and the system undergoes a solvable–unsolvable transition.

The above results allow us to determine how the integration step, $Δt$⁠, for the three integration schemes we have chosen, scales as the size of the problem increases. In Fig. 3, we define $Δtc$ such that $A(Δtc)=1/2$ and determine the relation between $Δtc$ and $N$⁠.

We find $Δtc$ scales as a power law with the number of variables $N$⁠. This is a major result because it shows that we do not need $Δt$ to decrease exponentially to have the same success rate. In Ref. 7, it was demonstrated, using topological field theory,^19,20 that in the ideal (noise-free) continuous-time dynamics, the physical time required for our dynamical system to reach the solution scales as $O(Nγ)$⁠, with $γ≤1$⁠.

Our numerical results show that, when simulating the dynamics with numerical integration schemes on classical computers, $Δtc∼N−δ$⁠. Coupled with the previous analytical results,⁷ this means that the number of integration steps scales as $O(Nγ+δ)$⁠. In other words, discretization only adds a $O(Nδ)$ overhead to the complexity of our algorithm, indicating that DMMs can be efficiently simulated on classical computers.²⁹ 

In Fig. 3, note that $Δtc>10−1$ for all three integration methods. This is quite unexpected for a stiff system simulated with explicit integration methods, because a large time step introduces large local truncation errors at each step of the integration. This error accumulates and should destroy the trajectory of the ODEs we are trying to simulate. However, we can still solve all tested instances with Eq. (2) even at such a large $Δt$⁠. In Sec. IV, we will provide an explanation of why this is possible.

Here, we note that choosing an appropriate $Δt$ could speed up the algorithm significantly, as smaller $Δt$ leads to excessive integration steps, while larger $Δt$ may cause the solver to fail. To find an appropriate time step for our scalability tests, we choose $Δt$ from Fig. 2 such that $95%$ of the trials have found a solution and multiply that $Δt$ by a “safety factor” of $0.6$ to make sure we hit the region of $A(Δt)≈1$⁠.

By employing the resulting value of $Δt$⁠, we plot in Fig. 4 the typical-case and 90th percentile scalability curves up to $N=105$⁠. For each variable size, $N$⁠, we perform 1000 3-SAT solution 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 Eq. (2) is evaluated], which differs from the number of integration steps by only a constant.

Both quantities show a power-law scalability for all three integration schemes. Since $Δtc$ is larger for higher-ordered ODE solvers, the latter ones typically require fewer steps to solve the problems. However, when taking the total number of function evaluations into account, the simplest forward Euler method becomes the most efficient integration scheme for our ODEs.

As anticipated, our results are unexpected in view of the large truncation errors introduced during the simulations of the ODEs of the DMMs. However, one might have predicted the numerical robustness upon considering the type of dynamics the DMM executes in phase space.

Instantons are topologically non-trivial solutions of the equations of motion connecting two different critical points (classical vacua in field-theory language).³⁰ In dissipative systems, as the ones considered in this work, instantons connect two critical points (such as saddle points, local minima or maxima) that differ in index (number of unstable directions). In fact, they always connect a critical point with higher index to another with lower index.^21,22

It was shown in Refs. 19 and 20 that, by stochastically quantizing the DMM’s equation of motion, one obtains a supersymmetric topological field theory, from which it can be deduced that the only “low-energy,” “long-wavelength” dynamics of DMMs is a collection of elementary instantons (a composite instanton). Critical points are also of topological character. For instance, their number cannot change unless we change the topology of phase space.³¹ 

In addition, given two distinct (in terms of indexes) critical points, there are several (a family of) trajectories (instantons) that may connect them, since the unstable manifold of the “initial” critical point may intersect the stable manifold of the “final” critical point at several points in the phase space. In the ideal (noise-free) continuous time dynamics, if the only equilibria (fixed points of the dynamics) are the solutions of the given problem (as shown, e.g., in Refs. 3 and 7), then the state trajectory is driven toward the equilibria by the voltages that set the input of the problem the DMM attempts to solve.

In the presence of (physical) noise instead, anti-instantons are generated in the system. These are (time-reversed) trajectories that connect critical points of increasing index. However, anti-instantons are gapped, in the sense that their amplitude is exponentially suppressed with respect to the corresponding amplitude of the instantons.³² This means that the “low-energy,” “long wave-length” dynamics of DMMs is still a succession of instantons, even in the presence of (moderate) noise.

Nevertheless, suppose an instanton connects two critical points, 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 original instanton has never been reached and that the critical point can be called absorbing, in the sense that the trajectory would “get stuck” on that one, or wanders in other regions of the phase space. In other words, in the presence of noise, there is a finite probability for some state trajectory to explore a much larger region of the phase space. The system then needs to explore some other (anti-)instantonic trajectory to reach the equilibria, solutions of the given problem. Nevertheless, due to the fact that anti-instantons are gapped, and the topological character of both critical points and instantons connecting them, if the physical noise level is not high enough to change the topology of the phase space, the dynamical system would still reach the equilibria.³³ 

On the other hand, if the noise level is too high, our system may experience a condensation of instantons and anti-instantons, which can break supersymmetry dynamically.³⁴ In turn, this supersymmetry breakdown can produce a noise-induced chaotic phase even if the flow vector field remains integrable. Although we do not have an analytical proof of this yet, we suspect that this may be related to the unsolvable phase we observe.

If we visualize the state trajectory as the one traced by a liquid in a corrugated landscape, the above suggests an intriguing analogy with the phenomenon of directed percolation (DP),³⁵ with the critical points acting as “pores” and instantons as “channels” connecting “neighboring” (in the sense of index difference) critical points.

DP is a well-studied model of a non-equilibrium (continuous) phase transition from a fluctuating permeable phase to an absorbing phase.³⁵ It can be intuitively understood as a liquid passing through a porous substance under the influence of a field (e.g., gravity). The field restricts the direction of the liquid’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 system goes through a phase transition from the absorbing phase into a permeable phase at a critical threshold $pc$⁠.

In numerical ODE solvers, discrete time also introduces some type of noise (truncation errors), and the considerations we have made above on the presence of anti-instantons still hold. However, numerical noise could be substantially more damaging than physical noise. The reason is twofold.

First, as we have already anticipated, numerical noise accumulates during the integration of the ODEs. In some sense, it is then always non-local. Physical noise, instead, is typically local (in space and time), hence it is a local perturbation of the dynamics.³⁶ As such, if it is not too large, it cannot change the phase space topology.

Second, unlike physical noise, integration schemes may change the topology of phase space explicitly. This is because, when one transforms the continuous-time ODEs (2) into their discrete version (a map), this transformation can introduce additional critical points in the phase space of the map, which were not present in the original phase space of the continuous dynamics. These extra (undesirable) critical points are sometimes called ghosts.²⁵ For instance, while the forward Euler method can never introduce such ghost critical points, irrespective of the size of $Δt$ (because both the continuous dynamics and the associated map have the same flow field $F$⁠), both the trapezoid and Runge-Kutta fourth order may do so if $Δt$ is large enough. These critical points would then further degrade the dynamics of the system.

With these preliminaries in mind, let us now try to quantify the analogy between the DMM dynamics in the presence of numerical noise and directed percolation.

First of all, the integration step, $Δt$⁠, must be inversely related to the percolation probability $p$⁠: when $Δt$ tends to zero, the discrete dynamics approach the ideal (noise-free) dynamics for which $p→1$⁠. In the opposite case, when $Δt$ increases, $p$ must decrease.

In directed percolation, the order parameter is the density of active sites.³⁵ This would correspond to the density of “achievable” critical points in DMMs. However, due to the vast dimensionality of their phase space (even for relatively small problem instances), this quantity is not directly accessible in DMMs.

Instead, what we can easily evaluate is the ratio of successful solution attempts: starting from a random point in the phase space [initial condition of the ODEs (2)], and letting the system evolve for a sufficiently long time, a path would either successfully reach the solution (corresponding to a permeable path in DP) or fail to converge (corresponding to an absorbing path in DP) . By repeating this numerical experiment for a large number of trials, the starting points of the permeable paths essentially fill the basin of attraction of our dynamical system. The advantage of considering this quantity is that the ratio of permeable paths in bond DP models can be calculated analytically.

Now we setup the correspondence of paths between DP and DMM. In DP, a permeable path is defined to be a path that starts from the top of the lattice and ends at the bottom and corresponds to a successful solution attempt in DMMs. Similarly, an absorbing path starts from the top and terminates within the lattice in DP and corresponds to a failed solution attempt in DMMs. Consider then a 3-SAT problem with only one solution. A DMM for such a problem can reach the solution from several initial conditions. This translates into a $D$-dimensional cone-shaped lattice for the DP model (see Fig. 5). The (top) base of the cone would represent the starting points, and the apex of the cone would represent the solution point. A permeable path connects the base to the apex, and an absorbing path ends in the middle of the lattice with all bonds beneath it disconnected.

Note that, in the DP model, it is possible to have different permeable paths starting from the same initial point. However, from the perspective of DMMs, there is no randomness in the dynamics, and each initial condition corresponds to a unique path. Since DP is a simplified theoretical model to describe numerical noise in DMMs, in this model we assume that the trajectory “randomly chooses” a direction when reaching a critical point. This is a reasonable assumption since there can be a large number of instantons starting from some critical point in the phase space, and it is almost impossible for us to actually monitor which instantonic path the trajectory follows.

With this correspondence in mind, the ratio of permeable paths in DP is analogous to the ratio of successful solution attempts in DMMs. Below, we will calculate this quantity exactly in DP and use this calculation to model the solvable–unsolvable transition we found in DMMs.

To begin with, we assume that all starting points are occupied with equal probability and calculate the expectation of the number of permeable and absorbing paths. This can be done iteratively: the expected number of permeable paths at a given site equals the sum of permeable paths at its predecessor sites, times $p$⁠, the percolation probability. Assuming the number of time steps required to reach the apex is $T$⁠, then the expected number of permeable paths

An illustration of the calculation on a cone-shaped 2D lattice is shown in Fig. 5.

The number of absorbing paths $⟨na,T⟩$ is trickier to compute. Details of the calculation can be found in  Appendix B, here we only give the approximate result (valid for $Dp>1$⁠) near the transition,

The ratio of permeable paths $r$ is simply $⟨np,T⟩/(⟨np,T⟩+⟨na,T⟩)$⁠,

We observe that the transition occurs at $pc∼eD$⁠. Near the transition, let us define $p=e+δD$⁠, where $δ$ is small. Then, $ln⁡Dp≈1+δe$⁠. Further using $limx→∞(1+1x)x=e$⁠, to order $O(δ)$⁠, Eq. (7) becomes

In the limit of $D→∞$⁠, the divergence comes from the $D$ in the exponent. Therefore, the transition happens exactly at $δ=0$⁠. When $δ>0$⁠, $r→1$⁠; when $δ<0$⁠, $r→0$⁠.

Note that when $δ≪0$⁠, the $erfc$ term dominates over the $e−Dδ/e$ term instead. However, in this case, $p$ is too small, and some approximations we made in  Appendix B to derive $⟨na,T⟩$ no longer hold. In this sense, Eqs. (7) and (8) are only valid near the transition point $pc=eD$⁠.

The behavior of Eq. (7) is plotted in Fig. 6 for different dimensions $D$⁠. Comparing to Fig. 2, we can already see a few similarities: both figures exhibit a sigmoidal behavior. However, as $Δt$ in DMMs is inversely related to $p$ in DP, they are curving in opposite directions.

Recall that the size of the basin of attraction, $A$⁠, in DMMs corresponds to the ratio $r$ in DP. To model their relation, we then use the ansatz $δ≡Dp−e=a(1NΔt−b)$⁠, with $a$ and $b$ being some real numbers, and fit $A$ to $Δt$ using Eq. (7). (Note that this trial function has a meaning only near $Δtc$⁠.)

We find that we can fit the curves reasonably well by fixing $a=5$⁠, and Fig. 7 shows the fitting result of $b$ and $D$⁠, where the number of variables $N$ ranges between $102$ and $104$⁠, and each data point is obtained by numerically integrating Eq. (2) for 1000 3-SAT solution attempts (100 instances and 10 initial conditions each), until a plateau, as in Fig. 1, has been reached. Note that the horizontal axis represents $1/(NΔt)$⁠.

The fitted parameters $b$ and $D$ are shown in Fig. 8, where both parameters exhibit a power-law scaling. At the transition threshold $δ=0$⁠, we have $b=1NΔt$⁠. Therefore, the power-law scaling of $b$ is closely related to the power-law scaling of $Δtc$ in Fig. 3, and these two scalings show similar trends.

The dimensionality $D$ is proportional to the dimensionality of the composite instanton (namely, the number of elementary instantons to reach the solution). In Ref. 20, this dimensionality was predicted to scale (sub-)linearly with the number of variables $N$ (in the noise-free case), and it is smaller than the dimensionality of the actual phase space [$(1+2×αr)N=17N$ in the present case]. This prediction agrees with our numerical results.

Finally, using the correspondence between DMMs in the presence of noise and DP, we can qualitatively explain why the DMMs numerically integrated still provide solution to the problem they are designed to solve even with such large numerical errors introduced during integration. In DMMs, the dimensionality of the composite instanton, $D$⁠, is usually very large, and DP tells us that the percolation threshold $pc=e/D∼e/N$⁠.

Therefore, even if most of the paths in the phase space are destroyed by noise, we can still achieve the solution as long as the probability of an instantonic path is larger than $e/D∼e/N$⁠. This argument, in tandem with the fact that critical points and instantons have a topological character,^21,22 ensures the robustness of DMMs in numerical simulations.

In conclusion, we have shown that the dynamics of digital memcomputing machines (DMMs) under discrete numerical solvers can be described as a directed percolation of state trajectory in the phase space. The inverse time step, $1/Δt$⁠, plays the role of percolation probability, and the system undergoes a “solvable-unsolvable phase transition” at a critical $Δtc$⁠, which scales as a power law with problem size. In other words, for the problem instances considered, we have numerically found that the integration time step does not need to decrease exponentially with the size of the problem, in order to control the numerical errors.

This result is quite remarkable considering the fact that we have only employed explicit methods of integration (forward Euler, trapezoid, and Runge–Kutta fourth order) for stiff ODEs. (In fact, the forward Euler method, although having the largest numerical error, solves the instances in the least amount of function evaluations.) It can be ultimately traced to the type of dynamics the DMMs perform during the solution search, which is a composite instanton in phase space. Since instantons, and the critical points they connect, are of topological character, perturbations to the actual trajectory in phase space do not have the same detrimental effect as the changes in topology of the phase space.

Since numerical noise is typically far worse than physical noise, these results further reinforce the notion that these machines, if built in hardware, would be topologically protected against moderate physical noise and perturbations.

However, let us note that we did not prove that the dynamics of DMMs with numerical (or physical) noise belong to the DP universality class. In fact, according to the DP-conjecture,^37,38 a given model belongs to such a universality class if (i) the model displays a continuous phase transition from a fluctuating active phase into a unique absorbing state; (ii) the transition is characterized by a non-negative one-component order parameter; (iii) the dynamic rules are short-ranged; and (iv) the system has no special attributes such as unconventional symmetries, conservation laws, or quenched randomness.

It is easy to verify that DMMs under numerical simulations satisfy properties (iii) and (iv). However, verifying properties (i) and (ii) is not trivial, as the relevant order parameters are not directly accessible due to the vast phase space of DMMs. Still, the similarities we have outlined in this paper, between DMMs in the presence of noise and DP, help us better understand how these dynamical systems with memory work and why their simulations are robust against the unavoidable numerical errors.

We thank Sean Bearden for providing the MATLAB code to solve the 3-SAT instances with DMMs, as well as insightful discussions. The work was supported by NSF under Grant No. 2034558. All memcomputing simulations reported in this paper have been done on a single core of an AMD EPYC server.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we show that the results presented in the main text hold also for other clause-to-variable ratios. As examples, we choose $αr=6$⁠. (Similar scalability results have been already reported in Ref. 7 for $αr=4.3$⁠.) In particular, we find similar scaling behavior for $Δtc$ for all clause-to-variable ratios, with of course, different power laws.

Figure 9 has been obtained as discussed in the main text and shows $Δtc$ vs number of variables for $αr=6$⁠.

In Fig. 10, instead we show the scalability curves for $αr=6$⁠, considering all three explicit integration methods. As in the main text, we find that the forward Euler method, although having the largest numerical error, solves the instances in the least amount of function evaluations for $αr=6$⁠. Every data point in the curves is obtained with 100 3-SAT instances, with 10 solution trials for each instance.

Here, we give a detailed calculation of the number of absorbing paths we outlined in Sec. IV E. We use $D$ to denote the dimension of the lattice where percolation takes place (corresponding to the dimension of the DMM phase space), $T$ to denote the number of steps to reach the bottom of the lattice (corresponding to the number of instantonic steps to reach the solution in DMMs), and $p$ to denote the percolation probability. Throughout the calculation, we use the approximation that $D→∞,T→∞$⁠, as both $D$ and $T$ are very large in DMMs.

As we illustrated in Sec. IV E, the expected number of permeable paths at time step $i$ (i.e., the $i$th level in the lattice in Fig. 5, where $i$ starts from 0) is

The number of absorbing paths at each site is the number of permeable paths at that site, times $(1−p)D$⁠, the probability that all bonds beneath it are disconnected. Note that this probability has a different expression at the boundary of the lattice, but as $D$ and $T$ are large, we ignore the boundary effect here. Meanwhile, the number of sites $mi$ at each time step $i$⁠, to the leading order, equals to the volume of a $(D−1)$-dimensional hyperpyramid,

Then, the total number of absorbing paths is

To get rid of the summation in Eq. (B3), we approximate it by adding the negligible $N=1$ term and replace the summation with an integral,

where $γ(s,x)=∫0xts−1e−tdt$ is the lower incomplete gamma function, and the result is valid for $Dp>1$⁠.

To further simplify $⟨na,T⟩$⁠, we use the relation³⁹ 

where $en(x)=∑k=0nxkk!$ is the truncated exponential series. Then,

where

is a number between 0 and 1.

Let us define

Using Stirling’s formula, we have

We can check that $g(k,x)$⁠, as a function of $k$⁠, reaches its maximum near $k=x$⁠. Let us expand $ln⁡k$ near $k=x$⁠,

Then,

Near the maximum $k=x$⁠, we have

which is a normal distribution with mean value $x$ and standard deviation $x$⁠. Figure 11 shows the comparison of the original $g(k,x)$ and its approximation. We can see that the approximation is very good.

Then, we have

where $Φ(x)=12π∫−∞xe−t2/2dt$ is the cumulative distribution function of the standard normal distribution.

Back to Eq. (B7), we finally have

where $erfc$ is the complementary error function.

In DMMs, each instanton connects two critical points whose indices (number of unstable directions) differ by 1.⁷ Since the index of a critical point equals at most the dimension $D$⁠, and the index of the equilibrium point is $0$⁠, we have $T+1=D$⁠, and Eq. (B14) simplifies to

Therefore,

which is Eq. (6) in the main text.