Computational phase transition signature in Gibbs sampling

Complexity theory aims to classify computational problems by their degree of difficulty—how much resource usage to solve the problems computationally. The satisfiability (SAT) problem was one of the first problems proven to be in a complexity class called NP-complete [1] and has been a useful workhorse for studying the properties of NP-complete problems. The SAT problem determines if a Boolean formula can evaluate to TRUE, and hence be satisfiable, or if it is otherwise unsatisfiable. A SAT formula presented in conjunctive normal form (CNF) consists of clauses connected by logical AND, with each clause being a combination (logical OR) of Boolean variables or their negations. When there are k variables in each clause, the problem is referred to as k-satisfiability (k-SAT). Considering random k-SAT where each clause is selected uniformly, SAT exhibits a rapid change from satisfiable instances to unsatisfiable instances at a certain clause to variable ratio (clause density) [2–11]. In addition, computer algorithms solving problem instances require increasing computational resources at this transition point. Such abrupt threshold behavior is known as the algorithmic or computational phase transition.

As information is necessarily represented in physical media, the processing, storage and manipulation of information is governed by the laws of physics. Indeed, the theory of computation is intertwined with the laws governing physical processes [12]. Many variants of naturally occurring processes (such as protein folding) have been shown to represent computationally significant problems, such as NP-hard optimization problems. Viewed another way, many physical systems and physical processes can be made to represent, and solve computational problem instances. In this regard, there are attempts to build contemporary physics based processors that solve problems instances by employing a physical process, such as annealing, quantum annealing and cooling [13–24].

Since the physical processes and the theory of computation are interwoven, we argue that signature of the algorithmic phase transition has potential to be observed in physical systems employing physical processes to solve problems instances. In particular, many physical systems have been build to solve computational problems by utilizing a cooling process which eventually leads a system down to equilibrium at lower temperature. In statistical mechanics, this equilibrium state is ideally described by the Gibbs state. In this regard, we consider the complexity in sampling and recovering the ground state from the Gibbs distribution of such physical systems.

Considering the Gibbs states for 2 (and 3)-SAT Hamiltonians, we calculated the ground-state occupancy as a function of clause density. We found that the algorithmic phase transition has a statistical signature in the Gibbs states of problem Hamiltonians generated randomly across the critical point. From our empirical observation, the ground-state occupancy of Gibbs states decrease for instances starting at the critical point of 2 (and 3)-SAT for fixed inverse temperatures and system sizes. This was confirmed by exact calculations of up to 26 binary units (spins) on a mid-scale supercomputer. Physical observation of the effect is hence within reach of near term and possibly even existing physical computing hardware: such as annealers [13], Ising machines [15, 16, 25] and quantum enhanced annealers [17, 18, 24, 26–30].

We also explore in this work a physical process such that the Gibbs distributions could be recovered in the long time limit. To simulate such physical systems, we consider the kinetic Ising model where the system couples with stochastic dynamics and is placed in contact with a heat bath [31–35]. A particular example of the kinetic Ising model that evolves by single-spin-flip dynamics is known as the Glauber model [31]. This model has been applied in the study critical behavior of Ising-like models [31, 36–39] and chaotic behavior [40]. In this setting, we simulated the time evolution which approximately represents the average behavior of spins with Glauber (spin–flip) stochastic dynamics [36], then numerically estimated the success probability of the simulated Glauber dynamics to sample the ground energy within a given error tolerance as a function of clause density. Considering random 2-SAT instances, we observed that the probability decreases for instances with density higher than the critical density.

The SAT problem is a decision problem to determine whether there exists a truth assignment to Boolean variables such that a given Boolean formula F evaluates to TRUE. If yes, we say F is satisfiable, otherwise unsatisfiable. The SAT problem was the first problem shown to be NP-complete by Cook in 1971 [1]. Many computational problems are known to be NP-complete through a polynomial mapping (Karp reduction) onto the k-SAT problem [41].

For example, considering a 3-SAT problem instance in CNF having two clauses over five variables,

$$\begin{align} \left(x_{1} \lor \neg x_{2} \lor x_{3}\right)\land \left(x_{1} \lor x_{4} \lor \neg x_{5}\right), \end{align} \tag{ 1 }$$

this instance is satisfiable as $x_1 = 1,\,x_2 = 1,\,x_3 = 0,\,x_4 = 0,\,x_5 = 0$.

Let F be a k-SAT instance consisting of M clauses in k CNF over N Boolean variables, the clause density of the instance is defined by the simple faction $\alpha = M/N$. For the above instance (1), $\alpha = 2/5$. We call F an instance of the random k-SAT problem if its clauses are randomly generated by uniformly selecting unique k-tuples from the union of a variable set (cardinality N > k) and its element wise negation.

The random k-SAT problem exhibits an abrupt threshold behavior separating a SAT phase from an UNSAT phase at the critical clause density $\alpha_c(k)$ [2–11, 42]. Moreover, it is observed that that the computational resources required for solving k-SAT instances start to increase from the transition point $\alpha_c(k)$.

For random 3-SAT instances, an abrupt transition in the probability of SAT has been numerically observed (though not formally proven) to be around clause density $\alpha_{c}(3) \approx 4.27$ [5–7, 10]. Most instances on the left (right) of the transition are satisfiable (unsatisfiable). It is statistically observed that computational time of decision SAT solvers peaks at the same point where hard instances for the 3-SAT problem are expected to concentrate (see figure 1 top right).

Zoom In Zoom Out Reset image size

The computational phase transition has better theoretical footing when considering the random 2-SAT problem, which can be solved in polynomial time [43, 44]. Both theoretical and numerical results [2, 8, 45] firmly position the transition at clause density $\alpha_{c}(2) = 1$ (see figure 1 top left). The further advantage is that additional slack bits are not required to embed 2-SAT into the Ising model. In other words, 2-SAT instances can be embedded directly into 2-body generalized Ising Hamiltonians (see appendix A).

We randomly generate 1000 instances of 2 (and 3)-SAT then calculate the percent of SAT (see figure 1 middle) and run-time (see figure 1 bottom) using a SAT solver, PycoSat package in Python 3 [46]. The sharp transition at the critical point can be observed in infinite limit of a ratio. In practice, finite sized effects broaden the transition region.

For 2-SAT, the theoretical results [45] establish a finite-size scaling window for the transition region. For $0 \lt \delta \lt 1$, this window is defined by $W(N,\delta) = (\alpha_{-}(N,\delta),\alpha_{+}(N,\delta))$, where $\alpha_{-}(N,\delta) = \sup(\alpha)$ such that the probability of satisfiablity is greater than $1-\delta$, similarly, $\alpha_{+}(N,\delta) = \inf(\alpha)$ such that the probability of satisfiablity is less than δ. As the width of the window tends to zero, $\alpha_{\pm}(N,\delta) \longrightarrow 1$, as $N \longrightarrow \infty$ for all δ.

For 3-SAT, it is still unknown where the exact transition point is. However, upper and lower bounds have been proven. The best lower bound at this time appears to be $\alpha_{lb} = 3.52$ [47] and the best upper bound is $\alpha_{ub} = 4.453$ [48].

In figure 1(middle and bottom), we can see that, as the number of variables increase, the transition sharpens and becomes closer to clause density 1 for 2-SAT (left) and 4.27 for 3-SAT (right).

Although the decision k-SAT problems exhibit the easy-hard-easy transition in the computational resources (see figure 1 bottom), the maximum k SAT (MAX-k-SAT) problem, the optimization variant of k-SAT where one seeks to determine variable assignments that maximize the number of satisfied clauses in a given instance, shows a different type of phase transitions, the easy-hard transitions [49, 50], which reflects the difference in computational complexity of these two problems. MAX-2-SAT exhibits a phase transition at the same critical clause density as 2-SAT [50].

We consider physical systems that encodes solutions of 2-SAT instances into the ground-state subspace of their corresponding Hamiltonians. This is done by using the standard Ising embedding [51–53] which embeds k-SAT instances directly into k-body generalized Ising Hamiltonians (see appendix A for details).

Considering the corresponding Gibbs state of 2-SAT Hamiltonians, we calculate the ground-state occupancy of such Hamiltonians from their Gibbs distributions as a function of clause density. This procedure can be formulated as follows.

Let $|i\rangle$ denote (possibly d-degenerate) lowest eigenstates of Hamiltonian $\mathcal{H} \unicode{x2A7E} 0$. We label these possibly degenerate states by letting i range from 1 up to d and call $E_\textrm{min}$ the lowest eigenvalue of $\mathcal{H}$

$$\begin{align} \forall i \in\left\{1, {\ldots}, d\right\}, ~\min_{\psi \in \mathbb{C}_2^N} \langle\psi|\mathcal{H}|\psi\rangle = \langle i|\mathcal{H}|i\rangle = E_\textrm{min} \unicode{x2A7E} 0. \end{align} \tag{ 2 }$$

In statistical mechanics, a physical system at thermal equilibrium is ideally described by a Gibbs state,

$$\begin{align} \rho_\beta = \frac{e^{-\beta \mathcal{H}}}{\mathcal{Z}}, \end{align} \tag{ 3 }$$

where β is inverse temperature and $\mathcal{Z}$ is the partition function defined as

$$\begin{align} \mathcal{Z} = \textrm{tr}\left\{e^{-\beta \mathcal{H}}\right\}. \end{align} \tag{ 4 }$$

Considering the Gibbs state (3), the occupancy in the (possibly d-degenerate) lowest-energy subspace of Hamiltonian $\mathcal{H}$ for fixed inverse temperature β becomes

$$\begin{align} P\left(\beta, N, \alpha \right) = \sum_{i = 1}^d\langle i|\frac{e^{-\beta\mathcal{H}}}{\mathcal{Z}}|i\rangle = \frac{d}{\mathcal{Z}}e^{-\beta E_\textrm{min}}. \end{align} \tag{ 5 }$$

The quantity $P\left(\beta, N, \alpha \right)$ in equation (5) provides our means to access the complexity of sampling the right solution from a thermal distribution. In the zero temperature limit, we have

$$\begin{align} \lim_{\beta \rightarrow \infty}P\left(\beta, N, \alpha \right) = \lim_{\beta \rightarrow \infty}\frac{1}{1+\frac{1}{d}\sum\limits_{E_j>E_\textrm{min}}e^{-\beta\left(E_j-E_\textrm{min}\right)}} = 1. \end{align} \tag{ 6 }$$

The average occupancy $P\left(\beta, N, \alpha \right)$ as a function of clause density α with different number of inverse temperature β is illustrated in figure 2(left). For each clause density, we randomly generated 1000 2-SAT instances from 26 variables. We then applied the embedding procedure from appendix A and calculated $P\left(\beta, N, \alpha \right)$.

Zoom In Zoom Out Reset image size

We observe that, on average, the ground state occupancy $P\left(\beta, N, \alpha \right)$ decrease for instances starting from the critical point ($\alpha \gt \alpha_c$). The same feature can be seen when varying the system size $N = \{5,10,15,20\}$ (see figure 3 left).

Zoom In Zoom Out Reset image size

The complexity landscape of $P\left(\beta, N, \alpha \right)$ is divided into two regions due to the density-dependent behavior of the ground subspace of SAT Hamiltonian. For $\alpha \lt \alpha_c$, the ground state energy is zero with its degeneracy decreases as clause density increases, which makes $P\left(\beta, N, \alpha \right)$ decrease with clause density. For $\alpha \gt \alpha_c$, the ground energy is an increasing function of α. However its degeneracy is almost constant with density α, hence causing $P\left(\beta, N, \alpha \right)$ to change its behavior around the critical density. Data for the density dependence of ground state energy and its degeneracy of 2-SAT are provided in appendix B, see figure B1.

One might be interested in the following question: on average, what inverse temperature is required to ensure that the occupancy of the ground thermal state, $P\left(\beta, N, \alpha \right)$, is greater than some fixed probability $\delta \in [0,1]$?

In figure 2(right), we plot the minimum value of inverse temperature ($\beta_\mathrm{min}$) satisfying the above condition with δ = 0.9 versus clause density for 2-SAT instances with 26 variables. On average, the instances starting from the critical point ($\alpha \gt \alpha_c$) require lower temperatures for sufficient ground state occupancy than ones with density $\alpha \lt \alpha_c$. The same feature also appears for other chosen values of $\delta = \{0.9, 0.7$, $0.5\}$ (see figure 3 right). For physical systems preparing in the Gibbs distributions, it informs us of the significance that temperature plays in obtaining above some percent ground-state occupancy.

Considering physical systems employing physical processes to solve problems instances, resultant distributions could deviate from the ideal Gibbs state. In this section, we study the complexity in sampling and recovering ground states in such physical systems. We simulated the Glauber time evolution [36] on the Ising spin models realizing 2-SAT instances, and numerically estimated the success probability in finding ground-states of 2-SAT Hamiltonians as a function of clause density.

The Glauber dynamics have been successfully used to treat and understand critical behavior in physics of complex systems [31, 36–40]. In the Glauber dynamics, the ith spin variable, $\sigma_i \in \{+1,-1\}$, is selected randomly at a time, and the spin flips its sign ($\sigma_i \rightarrow -\sigma_i$) at a transition rate that depends on the change in the energy of the system as a result of the flip.

For physical reasoning, ones assume that each transitioning process is in equilibrium with its reverse process. The equilibrium Boltzmann weight $e^{-\beta \mathcal{H}(\sigma)}/\mathcal{Z}$ is assigned with respect to the energy $\mathcal{H}(\sigma)$ of σ-configuration. According to this basic principle, we are able to fix the transition rates of the spin flipping by the the detailed balance condition:

$$\begin{align} \frac{w_i\left(\sigma \rightarrow \sigma^{^{^{\prime}}}\right)}{w_i\left(\sigma^{^{^{\prime}}} \rightarrow \sigma\right)} = \frac{e^{-\beta \mathcal{H}\left(\sigma^{^{^{\prime}}}\right)}} {e^{-\beta \mathcal{H}\left(\sigma\right)}}, \end{align} \tag{ 7 }$$

where $\sigma = \sigma_1 \cdots \sigma_i \cdots \sigma_N \in \{+1,-1\}^N$ denote a spin configuration of a system, and $\sigma^{^{^{\prime}}} = \sigma_1 \cdots (-\sigma_i) \cdots \sigma_N$ denote the configuration derived from σ due to flipping the ith spin: $\sigma_i \rightarrow -\sigma_i$.

From here onward, we consider the time evolution of Glauber dynamics [31, 36–38]. Let $\sigma_i \in \{+1,-1\}$ be the ith spin variables where $i = \{1, \cdots, N\}$. The Glauber time evolution is given by the following rate equation (8) which assigns the expected value of ith spin $\langle \sigma_i (t)\rangle \in [-1,+1]$ at time $t \in [0, \infty)$ as

$$\begin{equation} \frac{d}{dt} \langle \sigma_i \left(t\right)\rangle = -2 \langle \sigma_i \left(t\right) w_i\left(t\right)\rangle, \end{equation} \tag{ 8 }$$

with the transition rate

$$\begin{equation} w_i\left(t\right) = \frac{1}{2} \left(1-\tanh \frac{1}{2} \beta \delta \mathcal{H}_i \left(t\right) \right) \end{equation} \tag{ 9 }$$

satisfying the detailed balance condition (7) and $\delta \mathcal{H}_i (t)$ is the change of energy of a single-spin flip at the ith site.

Consider a 2-SAT instance having $M = \alpha N$ clauses, the corresponding Hamiltonian is expressed in the form of 2-body generalized Ising model

$$\begin{equation} \mathcal{H} = \frac{1}{4}\left[M \unicode{x1D7D9} + \sum_{i = 1}^{N} h_i \sigma_i +\sum_{j>i} J_{ij} \sigma_i \sigma_j\right] \end{equation} \tag{ 10 }$$

with appropriate coefficients $h_i, J_{ij}$ such that the ground state subspace is spanned by solutions of the problem instance, and the ground energy ($E_\mathrm{min}$) is equal to minimal number of unsatisfiable clauses. A detailed construction of such Hamiltonian is shown in appendix A. Thus the change of energy (10) by a single-spin flip at the ith site is given by

$$\begin{equation} \delta \mathcal{H}_i = -\frac{1}{2} \sigma_i \left( h_i + \sum_{j \neq i} J_{ij}\sigma_{j} \right), \end{equation} \tag{ 11 }$$

where $\sum_{j \neq i}$ is the summation over the jth spins that interact with the ith spin.

We will approximate the Glauber time evolution (8) by consideration of the limiting case where the equilibrium mean-field theory applies: the interaction is weak and long range. Substitute (11) into (9) and by using the fact that $\tanh(a \sigma_i) = \sigma_i \tanh(a)$ if $\sigma_i = \pm 1$, and $\sigma_i^2 = 1$, the Glauber rate equation (8) becomes the following approximate equation

$$\begin{equation} \frac{d}{dt} \langle \sigma_i \left(t\right)\rangle = -\langle \sigma_i \left(t\right)\rangle - \tanh \langle r_i \left(t\right)\rangle, \end{equation} \tag{ 12 }$$

where

$$\begin{equation} r_i\left(t\right) = \frac{1}{4} \beta \left( h_i + \sum_{j \neq i} J_{ij}\sigma_{j}\left(t\right) \right). \end{equation} \tag{ 13 }$$

In this limit, r_i is the sum of a large number of contributions from many different parts of the system, we assume that the fluctuations of these contributions about their mean values are independent, so that a law of large numbers applies, then the fluctuations of r_i about its mean values will be small, so that in the limit $\langle \tanh{r_i} \rangle = \tanh \langle r_i (t)\rangle$ [36].

One might ask the following question: considering the approximate Glauber rate equation (12), what is the probability of sampling energy $E \unicode{x2A7D} E_\mathrm{min}+\epsilon$ for a given error tolerance $\epsilon$ > 0? The probability of sampling energy $E \unicode{x2A7D} E_\mathrm{min}+\epsilon$ is defined as a simple fraction

$$\begin{equation} \mathrm{Prob}\left(E ~|~ E \unicode{x2A7D} E_\mathrm{min}+\epsilon\right) = \frac{n_{\epsilon}}{n_\mathrm{total}}, \end{equation} \tag{ 14 }$$

where $n_{\epsilon}$ is the number of samples such that $E \unicode{x2A7D} E_\mathrm{min}+\epsilon$, and $n_\mathrm{total}$ is the total number of samples.

We numerically simulated the approximate Glauber rate equation (12) and set the time-step to be 500 with the change of time $\Delta t = 0.05$ for every time evolution. For each 2-SAT instance, 300 initial states are sampling from the uniform distribution: $\langle \sigma_i (t = 0) \rangle \sim \textrm{Uni} \{-1,+1\}$ for all $i \in \{1,\cdots, N\}$. In order to obtain the statistics of the energy, we evaluate the average value of energy E over 300 instances for each clause density.

In figure 4(left), we shows the average values of energy error $(E-E_\mathrm{min})$ across phase transition for different inverse temperature. The average values of energy error peak around the critical point, suggesting a reasonable range of error tolerance $\epsilon$ in the estimation of $\mathrm{Prob}(E ~|~ E \unicode{x2A7D} E_\mathrm{min}+\epsilon)$ for any fixed inverse temperature β. For example, $\epsilon \unicode{x2A7E} 0.8$ for β = 2.

Zoom In Zoom Out Reset image size

Therefore, in figure 4(right), we illustrate the $\mathrm{Prob}(E ~|~ E \unicode{x2A7D} E_\mathrm{min}+\epsilon)$ across the phase transition for β = 2 and error tolerance $\epsilon \unicode{x2A7E} 0.8$. We observe a decreasing in $\mathrm{Prob}(E ~|~ E \unicode{x2A7D} E_\mathrm{min}+\epsilon)$ for instances with density $\alpha \gt \alpha_c$ of 2-SAT.

We predict that the algorithmic phase transition has a signature in thermal Gibbs states, which has implications for contemporary physical computing systems. The computational phase transition signature has strong theoretical footing when considering the 2-SAT problem—a problem which offers a direct embedding into 2-body generalized Ising Hamiltonians. For $k \unicode{x2A7E} 3$, the k-SAT problem on the other hand requires additional slack bits to create the k-body projectors to embed the problem instances.

We randomly generated 2-SAT instances and embed these instances in the generalized Ising Hamiltonian such that the ground-eigenspace is spanned by solutions to the problem instance. Considering the Gibbs states of 2-SAT Hamiltonians, we calculated the occupancy in ground-state subspace as a function of clause density. We found that the occupancy appear to decrease starting from the phase transition point for 2-SAT (see figure 2 left and figure 3 left). This finding indicates that the complexity in sampling solutions of 2-SAT instances from the thermal Gibbs states changes at the critical density. This feature is also seen in our numerical calculation for 3-SAT, see appendix C for more details.

Moreover one infers that temperature of the Gibbs state is the fleeting computational resource. Indeed, one argues that for any fixed inverse temperature β, there exists problem instances that would require a significant number of samplings to recover the ground-state.

We also consider physical systems employing physical processes to solve problem instances. We simulated the time evolution of corresponding physical systems of 2-SAT instances using the Glauber dynamics. We observed that the sampling success of finding the ground state at fixed inverse temperature becomes lower for instances starting at the critical point (see figure 4 right).

Our prediction connects the computational phase transition, which exhibits an importantly unsolvable feature in the theory of computational complexity, with physical processes that naturally reveal this feature in observable quantities. Since these 2-SAT instances can be directly embedded onto the generalize Ising Hamiltonian, and recently Ising devices have been built with increasing programability, the computational phase transition signature has the potential to be physically observed. Such an experiment provides a benchmark to locate where hard problem instances are most likely.

The large calculations were performed on the Skoltech's supercomputer, Zhores [54]. The code is written in C with nested OpenMP parallelism. The Hamiltonian is expanded in the shared memory of each node and we compute in parallel several clauses. Multiple nodes run the procedure concurrently with different random sequence to sample the statistics with MPI collectives. This implementation can use all the operational memory of each server. In figure 2, we illustrate results for N = 26; the program run on 8 nodes to collect statistics for around 13 h.

We also used the SAT solver called PycoSat in Python 3 [46] for the demonstration of the percent of SAT and run-time for 2 (and 3)-SAT in figure 1.

In simulating the approximate Glauber rate equation (12), we set the time-step to be 500 with the change of time $\Delta t = 0.05$ for every time evolution. For each instance, we uniformly random 300 initial states $\langle \sigma_i (t = 0) \rangle \sim \textrm{Uni} \{-1,+1\}$, for all $i \in \{1,\cdots, N\}$ in order to obtain the statistics of the energy, then evaluate average energy E over 300 instances of random 2-SAT for each clause density. In appendix D, we show an example of the simulation.

We thank the Skoltech-HPC team for supporting our large calculations and computational programming.

All data that support the findings of this study are included within the article (and any supplementary files).

The authors declare no competing interests.

The code for generating the data will be available on reasonable request.

All authors conceived and developed the theory and design of this study and verified the methods. J B proposed the study and supervised its execution. H P, V A, I Z developed and deployed the parallel code to collect numerical data. All authors contributed to interpret the results and wrote the manuscript.

Finding ground energy and ground state configurations of physical systems such as spin glasses is NP-hard [55]. k-SAT instances can be directly mapped onto the Hamiltonians in the form of generalized Ising model, see e.g. [51, 52]. The standard procedure is as follows.

We convert logical bits $\{0,1\}$ to spins

$$\begin{equation} \textrm{Logical}\,\ 0 \mapsto |0\rangle, \quad \textrm{Logical}\,\ 1 \mapsto |1\rangle. \end{equation} \tag{ A.1 }$$

The penalty Hamiltonian $\mathcal{H}_\mathrm{SAT}$ is constructed by real linear extension of $\{P_0, P_1, \unicode{x1D7D9}\}$ by means of the following mapping between Boolean variables and projectors,

$$\begin{equation} x_j \longrightarrow P_j^0, \qquad \neg x_j \longrightarrow P_j^1, \end{equation} \tag{ A.2 }$$

and

$$\begin{equation} \wedge \longrightarrow +, \qquad \vee \longrightarrow \otimes, \end{equation} \tag{ A.3 }$$

where $P_j^{1} = |1 \rangle\langle 1|$ and $P_j^{0} = |0 \rangle\langle 0|$ acting on the jth spin.

By this construction, clauses in a SAT instance are mapped onto projectors $P_{i j \cdots k}^{\alpha \beta \cdots \gamma} = P_{i}^{\alpha} \otimes P_{j}^{\beta} \otimes \cdots \otimes P_{k}^{\gamma}$. Hence, the Hamiltonian $\mathcal{H}_\mathrm{SAT}$ can be constructed by summing over all the clauses in an instance,

$$\begin{equation} \mathcal{H}_\mathrm{SAT} = \sum_{l = 1}^{M} \mathcal{C}_{l} \lbrace P_{i j{\ldots}k}^{\alpha \beta{\ldots}\gamma} \rbrace, \end{equation} \tag{ A.4 }$$

where $\mathcal{C}_{l}$ assigns the value of $\alpha, \beta, \cdots, \gamma \in \{0, 1\}$ and the value of spin indices $i, j, \cdots, k \in \{1, 2, \cdots, N\}$ corresponding to the lth clause in the instance. The ground state space is spanned by solutions of the SAT problem, and the ground energy is equal to minimal number of unsatisfiable clauses.

Substituting the projectors $P^{\alpha}_j = \frac{1}{2}(\unicode{x1D7D9} + (-1)^{\alpha} \sigma^{z}_{j})$ in equation (A.4), where σ^z is the Pauli Z matrix, the k-SAT Hamiltonian can be written in the form of the generalized Ising Hamiltonian with at most k-body interaction.

In case of 2-SAT, the Hamiltonian $\mathcal{H}_{2SAT}$, which is a sum over projectors onto 2-SAT clauses, can be written in the form of 2-body Ising spin Hamiltonian,

$$\begin{align} \mathcal{H}_\mathrm{2SAT} & = \sum_{l = 1}^{M} \mathcal{C}_l \lbrace P^{\alpha}_j \otimes P^{\beta}_k\rbrace \nonumber\\ & = \frac{1}{4} \left( M \unicode{x1D7D9} + \sum_{j}h_{j} \sigma^{z}_{j} + \sum_{j \lt k}J_{jk}\sigma^{z}_{j}\sigma^{z}_{k}\right). \end{align} \tag{ A.5 }$$

The coefficient h_j indicates the local field of the jth spin and the couplings $\mathcal{J}_{jk}$ encode the 2-body interaction between the jth spin and the kth spin.

This method provides a way to physically realize SAT instances as spin Hamiltonians. If such a physical system is built, cooling the system into its ground state is equivalent to solving the problems.

In the following figure B1, the average values of (left) ground state energy and (right) its degeneracy over 300 random 2-SAT Hamiltonians of ten variables for clause density $\alpha \unicode{x2A7D} 6.$

Zoom In Zoom Out Reset image size

We present here the supplementary data for the random 3-SAT problems. In figure C1, we illustrate the occupancy $P\left(\beta, N, \alpha \right)$ for various inverse temperature (left) and the minimum inverse temperature $\beta_\mathrm{min}$ such that $P\left(\beta, N, \alpha \right) \unicode{x2A7E} 0.9$ (right) across the algorithmic phase transition for 26 spin variables.

Zoom In Zoom Out Reset image size

In figure C2, we show the occupancy for different problem size $N = \{5,10,15,20\}$ (left) and the minimum inverse temperature for various values of $\delta = \{0.9, 0.7$, $0.5\}$ across the algorithmic phase transition for N = 10 spins (right).

Zoom In Zoom Out Reset image size

Considering a 2-SAT instance with ten variables and ten clauses (clause density α = 1)

$$\begin{equation} \begin{split} &\left(\neg x_{6} \lor \neg x_{9}\right) \land \left(\neg x_{1} \lor x_{8}\right) \land \left(x_{7} \lor x_{10}\right) \land \left(x_{4} \lor \neg x_{9}\right) \land \left(x_{5} \lor x_{6}\right) \land\\ &\left(\neg x_{3} \lor x_{4}\right) \land \left(\neg x_{6} \lor x_{7}\right) \land \left(\neg x_{2} \lor \neg x_{6}\right) \land \left(x_{3} \lor \neg x_{9}\right) \land \left(\neg x_{3} \lor x_{6}\right) \end{split}. \end{equation} \tag{ D.1 }$$

We show the numerical convergence in expected spin variables and its corresponding energy in the simulation of the Glauber time evolution for different values of $\beta = \{1,2,3\}$ in figure D1.

Zoom In Zoom Out Reset image size