Computationally **challenging optimisation problems** have always been of special interest by various researchers around the globe. This is primarily due to them often having a very high dimensional search space, or having **highly complex **and non-linear objective functions at their core, which classical gradient-based methods fail to tackle efficiently. This has been the main reason for the development of **metaheuristic algorithms** that take inspiration from our surroundings (nature, swarms and physical processes) and can provide a computationally cheap yet robust optimisation procedure for such hard problems at hand. Parallely, researchers have also noticed the undeniable success of modeling physics’ processes to study highly complex phenomena, both in real-world and computer science. For instance, resource allocation problem has been well tackled by statistical mechanics models, while certain aspects of statistical thermodynamics have been employed to explain micro-evolution of species, and so on. As a result, in spite of swarm-inspired algorithms being in the forefront as robust optimisers, researchers have shown keen interest in adapting principles and theories of physics and applying them to solve real-world optimisation problems.

Recently, the world of metaheuristics has seen the advent of several **novel search mechanisms** based on various non-linear physics processes. The novelty of these approaches lies in the fact that the non-linear physical phenomena are leveraged as backbones to be modeled upon in order to formulate efficient search algorithms, whose mechanism is quite different from the conventional swarm and evolutionary algorithms. Such physics-inspired algorithms have shown great promise and robustness as **global optimisation strategies**.

A really compelling overview of the different approaches has been published here:

A Brief Overview of Physics-inspired Metaheuristic Optimization Techniques

Soumitri Chattopadhyay, Aritra Marik, Rishav Pramanik (2022)

**Digital Memcomputing Machines**

One of the approaches we are discussing today are Digital Memcomputing Machines. In “Memcomputing NP-complete problems in polynomial time using polynomial resources and collective states“, Memcomputing is introduced as a novel **non-Turing paradigm** of computation that uses interacting memory cells (memprocessors for short) to store and process information on the same physical platform. It is claimed that it was recently proven mathematically that memcomputing machines have the same computational power of **nondeterministic Turing machines**. Therefore, they should solve NP-complete problems in **polynomial time** and, using the appropriate architecture, with **resources** that only grow **polynomially** with the input size. Our implementation here reproduces the results as shown in the paper:

As can be seen, for the problem sizes tested, the authors find a power-law scaling for the median number of integration steps for the simulations of DMMs. They also find that integration time variable (t), CPU time, and long-term memory (xl) are bounded by a polynomial scaling, and the average step size shows power-law decay. Te optimised WalkSAT algorithm they have used instead exhibits an exponential scaling at relatively small problem sizes, confirming the previous results. An exponential scaling is also observed for the SID algorithm.

The Digital Memcomputing Machine is implemented as a simulation of an analog circuit, specifically by numerically integrating ODEs (Ordinary Differential Equations) on a traditional computer. At first, the CNF formulation of the SAT problem is constraint with the following equation:

Then, the equations of motion for the voltages and auxiliary variables to represent the Memristor based reversible gates are defined as followed:

We are implementing the solver in C++ using the Boost library, especially because of their provided ODEInt library which allows a fairly straight forward approach to implementing numerical integrations of ODEs.

The main steps of the solver are based on the following building blocks:

[1] Reading the CNF formulation

[2] Initialising the parameters and initial conditions

[3] Numerical integration of the ODEs

We assume that the reader is already familiar with SAT solving and CNF formulations, so we will start with the parameter initialisation:

**Initialising the parameters and initial conditions**

The equations of motion highly depend on a few parameters which have been defined in the paper as followed. We will find that changing (and tuning) of these parameters will have significant impact on the solver’s capability to reach the solution:

// growth rate for long term memory Xl TFloat dmm_alpha = TFloat(5.0); // growth rate for short term memory Xs dmm_beta = TFloat(20.0); // restriction for Cm in short term memory Xs TFloat dmm_gamma = TFloat(TFloat(25)/TFloat(100)); // restriction for Cm in long term memory Xl TFloat dmm_delta = TFloat(TFloat(5)/TFloat(100)); // remove spurious solution X,s,m = 0 TFloat dmm_epsilon = TFloat(TFloat(1)/TFloat(10)); // reduction factor of rigidity G (learning rate) TFloat dmm_zeta = TFloat(TFloat(1)/TFloat(10));

**Numerical integration of the ODEs**

We can see that the numerical integration of the ODEs can be realised by a loop going through each clause to calculate the derivate (dxdt) for the next time step:

// Loop ------------------------------------------------------------ //set all V to 0.00: fill( dxdt.begin() , dxdt.begin()+n , TFloat(0.0) ); // loop through each clause: for (int clause = 0; clause < m; clause++) { // number of literals in this clause: int ksat = clauseSizes[clause]; // Xl & Xs: TFloat Xs = x[clause+n]; if (Xs<0.0) Xs = TFloat(0.0); if (Xs>1.0) Xs = TFloat(1.0); TFloat Xl = x[clause+n+m]; if (Xl<1.0) Xl = TFloat(1.0); if (Xl>xl_max) Xl = TFloat(xl_max); TFloat C = TFloat(0.0); TFloat Ri, Rj, Rk, Gi, Gj, Gk; // 3-sat: if (ksat==3) { // +1 if literal is >0, otherwise -1 int Qi = (cls[clause*MAX_LITS_SYSTEM+0]>0)? 1:-1; // +1 if literal is >0, otherwise -1 int Qj = (cls[clause*MAX_LITS_SYSTEM+1]>0)? 1:-1; // +1 if literal is >0, otherwise -1 int Qk = (cls[clause*MAX_LITS_SYSTEM+2]>0)? 1:-1; int liti = abs(cls[clause*MAX_LITS_SYSTEM+0]); int litj = abs(cls[clause*MAX_LITS_SYSTEM+1]); int litk = abs(cls[clause*MAX_LITS_SYSTEM+2]); TFloat Vi = x[liti-1]; if (Vi<-1.0) Vi = TFloat(-1.0); if (Vi>1.0) Vi = TFloat(1.0); TFloat Vj = x[litj-1]; if (Vj<-1.0) Vj = TFloat(-1.0); if (Vj>1.0) Vj = TFloat(1.0); TFloat Vk = x[litk-1]; if (Vk<-1.0) Vk = TFloat(-1.0); if (Vk>1.0) Vk = TFloat(1.0); TFloat i = TFloat(1.0)-TFloat(Qi)*Vi; TFloat j = TFloat(1.0)-TFloat(Qj)*Vj; TFloat k = TFloat(1.0)-TFloat(Qk)*Vk; C = TFloat(fmin(i, fmin(j, k))); C = C / TFloat(2.0); if (C<0.0) C=TFloat(0.0); if (C>1.0) C=TFloat(1.0); // equation Gn,m(vn,vj,vk)= ... Gi = TFloat(Qi) * fmin(j,k) / TFloat(2.0); Gj = TFloat(Qj) * fmin(i,k) / TFloat(2.0); Gk = TFloat(Qk) * fmin(i,j) / TFloat(2.0); // equation Rn,m (vn , vj , vk ) = .... if (C != (TFloat(1.0)-TFloat(Qi)*Vi)/TFloat(2.0) ) {Ri = TFloat(0.0);} else {Ri = (TFloat(Qi) - Vi) / TFloat(2.0);} if (C != (TFloat(1.0)-TFloat(Qj)*Vj)/TFloat(2.0) ) {Rj = TFloat(0.0);} else {Rj = (TFloat(Qj) - Vj) / TFloat(2.0);} if (C != (TFloat(1.0)-TFloat(Qk)*Vk)/TFloat(2.0) ) {Rk = TFloat(0.0);} else {Rk = (TFloat(Qk) - Vk) / TFloat(2.0);} // equation Vn = ... TFloat _Vi = Xl * Xs * Gi + (TFloat(1.0) + m_zeta * Xl) * (TFloat(1.0) - Xs) * Ri ; TFloat _Vj = Xl * Xs * Gj + (TFloat(1.0) + m_zeta * Xl) * (TFloat(1.0) - Xs) * Rj ; TFloat _Vk = Xl * Xs * Gk + (TFloat(1.0) + m_zeta * Xl) * (TFloat(1.0) - Xs) * Rk ; // update voltages: dxdt[liti-1] = dxdt[liti-1] + _Vi; dxdt[litj-1] = dxdt[litj-1] + _Vj; dxdt[litk-1] = dxdt[litk-1] + _Vk; //update loc: if (C<0.5) loc--; //this clause is sat, reduce loc // Calculate new Xs: dxdt[n+clause] = m_beta * (Xs + m_epsilon) * (C - m_gamma); // Calculate new Xl: dxdt[n+m+clause] = dmm_alpha_cm[m_nodeid*m+clause] * (C - m_delta); } }

One of the potential pitfalls of the implementation is the fact that the variables used in the ODE are bounded. For example, the memory variables are bounded, with xs,m ∈ [0, 1] and xl,m ∈ [1, 104M]. Te boundedness of voltage and memory variables implies that there are no diverging terms in the equations. This can be seen in the code snipped above.

Using the Boost library, the integration itself can be done quite easily:

integrate_n_steps(stepper_euler, odeS(node->id, thread_params[node->id*PARAM_COUNT+0], thread_params[node->id*PARAM_COUNT+1], thread_params[node->id*PARAM_COUNT+2], thread_params[node->id*PARAM_COUNT+3], thread_params[node->id*PARAM_COUNT+4], thread_params[node->id*PARAM_COUNT+5]), x, t, thread_params[node->id*PARAM_COUNT+8], TFloat(maxsteps), write_ode(node->id));

This example uses a simple forward Euler integration scheme (which has been used by the authors of the paper), but we have also integrated Runge-Kutta and others.

### Results

To test the performance of our implementation, we are applying the same CNF instances the authors have been using in their publication. We start with a very simple instance, with N=100 and a ratio of 4.3. We run the solver with the following command:

./dmm -i transformed_barthel_n_100_r_4.300_p0_0.080_instance_001.cnf

The satisfying assignment is found in 0.01s. We gradually increase the complexity (N) to see if the solver actually scales linear with the problem size (and not exponentially). Here are the results (running on an MacBook 2021):

N=100 0.01s

N=1000 0.03s

N=10000 0.08s

N=100000 2.23s

Given that the processing time for each time step increases as the loop has to go through more and more clauses to calculate the time derivate, solutions to the SAT problems can be found effortless and instantly, ranging from N=100 to N=100,000, **confirming the results** of the paper.

We have presented an implementation of an **efficient dynamical-system** approach to solving Boolean satisfiability problems. Along with arguments for polynomial-time scalability in continuous time, we have found that the numerical integration of the corresponding ODEs show power-law scalability for typical cases of 3-SAT instances which required exponential time to solve with successful algorithms. Te efficiency derives from collective updates to the variables during the solution search (long-range order). In contrast to previous work, our dynamical systems do not suffer from exponential fluctuations in the energy function due to chaotic behaviour. Te dynamical systems we propose find the solution of a given problem without ever entering a chaotic regime, by virtue of the variables being bounded. The implication is that a hardware implementation of DMMs would only require a **polynomially-growing energy cost**. The work then also serves as a counterexample to the claim that chaotic behaviour is necessary for the solution search of hard optimisation problems. In fact, the authors find chaos to be an undesirable feature for a scalable approach. Although these analytical and numerical results do not settle the famous P vs. NP question, they show that appropriately designed physical systems are very useful tools for new avenues of research in constraint satisfaction problems.

The complete **source-code** of the implementation is available on our GitHub repository. More of our research and work in physics and bio inspired computing is also available here.