# Relaxation in graph coloring and satisfiability problems

## I Introduction

Computers have made it possible for physicists to do experiments without leaving their offices. By simulating systems from Nature, or simple models from theoretical physics, a tremendous amount of information can be gained. The advance of computer hardware and software has made it possible to simulate systems consisting of vast numbers of interacting particles. Physics, in turn, has inspired new algorithms for solving problems in computer science. For instance, the method of simulated annealing , based on Monte Carlo simulations with the Metropolis algorithm , can sometimes find good solutions much quicker than traditional algorithms .

Physicists have also started to study systems that are not found in Nature, but instead come from computer science . As an example, phase transitions in optimization problems have been discovered and studied using statistical mechanics . Other problems studied using physical methods include the knapsack-problem , graph-partitioning , minimax-games , the 8-Queens problem , number partitioning  and the stable marriage problem . Field theory has also been used to study, e.g., the enumeration of Hamiltonian cycles on graphs  and coloring of random, planar graphs .

Here we present a study of the $`T=0`$ relaxation behavior of hard optimization problems. Using Monte Carlo simulation, we have measured the energy of a system starting in a random (excited) state and slowly relaxing into the ground state or a low-lying excited state. We find qualitatively different behavior for hard and easy instances of the problems. One of many models (e.g., ) for which relaxation has been studied is the ferromagnetic Ising model on a regular lattice. The models studied here differ from the Ising model in several respects. They are random, i.e., there is no pattern in the interactions between the spins, and the interactions have infinite range. K-COL can be viewed as a Potts model on a random graph with finite connectivity. Another important difference is that these models can suffer from frustration.

This article is organized in the following way: Section II introduces some concepts from computer science, while the problems we study are described in section III. Previous work on the solvability transition is reviewed in section IV and some approximate explanations in section V. Section VI describes our results on the relaxation behavior, section VII compares the results to other local search methods. Additional measurements, e.g., of the fraction of persistent spins, are found in section VIII. Section IX discusses possible explanations for the relaxation behavior and comments on the importance of entropy barriers. Conclusions and a discussion are contained in section X.

## II Computer Science for Physicists

Computer scientists classify problems according to the maximal amount of resources needed for their solution. The most important resource is time, but it is also possible to distinguish between problems that require qualitatively different amounts of memory. For example, a list of $`N`$ elements can always be sorted in time less than $`kN\mathrm{log}N`$, where $`k`$ is some constant . The problems whose running time on a universal Turing machine (e.g., ) is bounded by a polynomial in their size are said to be in the class P. The important class NP (for non-deterministic polynomial) consists of those problems where it can be checked in polynomial time whether a proposed solution actually solves the problem. (A non-deterministic Turing machine would be able to solve NP problems in polynomial time.) It is obvious that $`𝐏\mathrm{𝐍𝐏}`$, but there is no proof that $`𝐏\mathrm{𝐍𝐏}`$. However, most people believe that there are NP problems whose worst-case instances take exponential time to solve on a universal Turing machine.

The class NP-complete (or NPC) are the most important problems in NP. A problem of size $`N`$ is in NPC if all other NP problems can be transformed into it in time at most polynomial in $`N`$. A method to solve an NPC problem efficiently can thus be used to solve any NP problem efficiently. It is known that if $`𝐏\mathrm{𝐍𝐏}`$ then there are problems in NP that are in neither P nor NPC. A problem is called NP-hard if it is at least as difficult as the most difficult NP problems; NPC is the intersection of NP and NP-hard.

It is worth emphasizing that it is the worst case complexity that determines whether a problem is in NPC. The average time needed to solve an NPC problem (given a distribution of problems) may still be polynomial in problem size. The properties of K-COL and K-SAT studied here are related to average case behavior, and do not address the question of whether $`𝐏\mathrm{𝐍𝐏}`$.

A nice introduction to most elementary concepts from theoretical computer science can be found in . A modern reference on complexity theory and NP problems is , while  has an extensive list of NPC problems.

## III K-SAT and K-COL

Two important problems in NPC are graph coloring (K-COL) and satisfiability testing (SAT). Graph coloring is the problem of coloring a graph with $`N`$ vertices and $`M`$ edges using $`K`$ colors so that no two adjacent vertices have the same color. In physical terms, K-COL is the problem of finding a ground state without frustrated bonds in an antiferromagnetic $`K`$-state Potts model on a random graph. Related models have been studied by Baillie, Johnston, and coworkers (e.g. ), who considered Potts models on $`\varphi ^n`$-model Feynman diagrams. They found similarities between models on $`\varphi ^3`$ and $`\varphi ^4`$ graphs and Bethe lattices, and showed that mean-field theories work well for describing both ferromagnetic and antiferromagnetic models on Feynman diagrams.

The most natural application of graph coloring is in scheduling problems. For example, a school where each teacher and student can be involved in several different classes must schedule the classes so that no collisions occur. If there are $`K`$ different time slots available, this is K-COL.

Satisfiability was the first problem shown to be in NPC . It is the problem of finding an assignment of true or false to $`N`$ variables so that a boolean formula in them is satisfied. In K-SAT, this formula is written in conjunctive normal form (CNF), that is, it consists of the logical AND of $`M`$ clauses, each clause being the OR of $`K`$ (possibly negated) variables, where the same clause may appear more than once in a formula. For example, $`(xy)(y\neg z)`$ is an instance of 2-SAT with two clauses and three variables. Applications of K-SAT include theorem proving, VLSI design, and learning.

In K-SAT, each clause forbids one of the $`2^K`$ possible assignments for its variables. In the same way, an edge in a graph forbids $`K`$ of the $`K^2`$ different colorings of its vertices. For both problems, there are $`M`$ constraints on the solutions. The energy $`ϵ`$ of a problem instance is defined as the number of unsatisfied constraints per variable. Both K-SAT and K-COL are in P for $`K=2`$ and in NPC for $`K3`$ . The related problem (MAX-K-SAT) of trying to minimize the number of unsatisfied clauses in K-SAT is in NPC even for $`K=2`$.

The most important problems are those where the number of constraints is of the same order as the number of variables, $`M=\alpha N`$. The scheduling problem described above fulfills this condition, for example. For graph coloring, $`\alpha =\gamma /2`$, where $`\gamma `$ is the connectivity of the graph. The connectivity (or average degree) is defined as the mean number of edges exiting each node. For a graph with $`n`$ vertices and $`e`$ edges, it is $`2e/n`$.

Below, we will use $`\alpha `$ for K-SAT and $`\gamma `$ when we talk of K-COL. We will concentrate on K-COL; most results for K-SAT are similar.

## IV The transition

For physicists, an interesting property of these problems is that they contain phase transitions . As the number of constraints increases, there is a transition from a region where almost all instances of the problem are solvable to a region where practically none can be solved. Physically, the transition is from a region where the ground state has zero energy to one where it is finite. For $`K=2`$ the transition can be seen as a transition between a P problem (finding a perfect solution to the problem) and an NPC problem (finding an assignment of variables that minimizes the number of unsatisfied constraints).

An approximation similar to mean-field theory has been used by Williams and Hogg (e.g., ) and others to explain some of the properties of the phase transition. Recently, Friedgut has also made some progress towards showing rigorously the existence of a sharp transition in solvability for both K-SAT  and K-COL .

Related to this phase transition in problem solvability, there is a transition in how difficult it is to solve a problem or show that no solutions exist . If there are few constraints on the solution (the problem is underconstrained), it is easy to find one. Similarly, if there are so many constraints that the problem is overconstrained, not much effort is needed to show that it is unsolvable. In between these regions, where the problems are critically constrained, there is a peak in problem difficulty. This is called the “easy-hard-easy” transition .

K-SAT has recently been studied by a number of physicists. Kirkpatrick and Selman  studied the phase transition using finite-size scaling methods, and Monasson and Zecchina  used the replica method  to show that the entropy of K-SAT stays finite at the transition. This means that below the transition there are several solutions to each problem, all of which develop inconsistencies as the critical $`\alpha `$ is passed. Another problem that has been studied using statistical mechanics is the number partitioning problem. Mertens has recently shown that it too has a phase transition. The relevant parameter here is the ratio between the number of bits of input data and the number of variables.

Another problem in NPC that also shows a transition  is the traveling salesperson problem (TSP), where the objective is to find a tour of minimum length visiting $`N`$ given distinct cities. A difficulty in studying this problem is that there is no natural parameter (like $`\alpha `$ and $`\gamma `$) that distinguishes between under- and overconstrained problems. To get one, the TSP must be reformulated as a decision problem: is there a Hamiltonian path of length less than $`l`$? The parameter $`l`$ plays the same rôle as $`\alpha `$ — for a given distribution of problems there is an $`l_c`$ such that if $`ll_c`$, almost all instances have a tour with length $`<l`$, but if $`ll_c`$ practically no such tours exist. Traditionally, most NPC problems are formulated as decision rather than optimization problems.

There are many NPC problems that contain no obvious parameter which makes it difficult to say if the solvability phase transition exists in all NPC-problems or in just a few. There have been attempts to formulate a more general parameter (e.g., ), with the drawback that it requires us to approximate the number of solutions.

Phase transitions have also been found in problems beyond NPC, e.g., in QSAT , a harder version of satisfiability where the boolean variables are quantified by either $``$ or $``$ (in ordinary SAT, all variables are existentially quantified). This problem is known to be PSPACE-complete , meaning that it is at least as hard as all problems that can be solved by a universal Turing machine without time limits but using memory at most polynomial in problem size.

## V An approximate theory for the transition

The approximation proposed by Williams and Hogg  assumes that the constraints in the problem are independent. In physical terms it simply means performing an annealed rather than a quenched average over the disorder. It is exact for graphs without loops and for satisfiability problems where no variable is contained in more than one clause. The probability that an independent constraint is violated is $`p=1/2^K`$ for K-SAT and $`p=K/K^2=1/K`$ for K-COL. The probability of having none of $`M`$ constraints violated can be approximated as $`(1p)^M`$, ignoring correlations between constraints, such as triangles in graphs. The number of solutions for K-COL is then

$$N_{\text{sol}}=K^N(1\frac{1}{K})^{\gamma N/2},$$
(1)   
and for K-SAT the expression is

$$N_{\text{sol}}=2^N(1\frac{1}{K^2})^{\alpha N}.$$
(2)   
Using the inclusion-exclusion principle it is possible to write an exact expression for $`N_{\text{sol}}`$ . The inclusion-exclusion principle is the generalization of the simple formula

$`P(AB)=P(A)+P(B)P(AB)`$

from mathematical statistics. If we let $`A_i`$ be the event that constraint $`i`$ is violated, it expresses the probability that any (i.e., at least one) constraint is violated in terms of the probabilities of one, two, three or more constraints being violated simultaneously

$$P(_iA_i)=\underset{r=1}{\overset{M}{}}(1)^{r+1}S_r,$$
(3)   
where $`S_r`$ is the probability of exactly $`r`$ constraints being violated simultaneously. The number of solutions can now be found as

$$N_{\text{sol}}=N_{\text{tot}}(1P(_iA_i)),$$
(4)   
where $`N_{\text{tot}}`$ is the number of possible assignments of the variables, $`N_{\text{tot}}=K^N`$ for K-COL and $`N_{\text{tot}}=2^N`$ for K-SAT. For K-COL, $`S_1=MK^1`$, since there are $`M`$ edges and each of them eliminates $`K^{N1}`$ (of the $`K^N`$) solutions. For $`S_2`$, we need to express the number of states that are eliminated by each of two edges. This is given by $`\left(\begin{array}{c}M\\ 2\end{array}\right)K^2`$, while the expression for

$$S_3=\left(\begin{array}{c}M\\ 3\end{array}\right)K^3+(K^2K^3)t,$$
(5)   
requires knowledge of the number of triangles, $`t`$, in the graph. Expression (5) can be understood by noting that if two edges in a triangle are frustrated, the third is always frustrated too. It can be shown that $`t`$ is Poisson-distributed with mean $`\gamma ^3/6`$. To calculate $`S_i`$ for $`i4`$, we also need to know the distribution of more complex sub-graphs.

The critical value of the parameter can be approximated as that $`\gamma `$ which gives $`N_{\text{sol}}=1`$ in (1), giving

$$\gamma _c=2\frac{\mathrm{log}K}{\mathrm{log}(1\frac{1}{K})}.$$
(6)   
For $`K=3`$, equation (6) gives $`\gamma _c=5.4`$ for K-COL and $`\alpha _c=5.2`$ for K-SAT. These values are larger than the experimental values of $`\gamma _c=4.6`$ and $`\alpha _c=4.17`$. For K-SAT, this approximation has been independently introduced several times \[39, and references therein\].

This calculation of the critical value of $`\gamma `$ ignores all correlations between different constraints in the problem. It gives an upper bound for $`\gamma _c`$ and is analogous to studying a forest, a graph without cycles, in which all edges are violated with a probability $`p`$. Taking correlations into account reduces the number of solutions.

Kirousis et al  have introduced a new method to get an upper bound for the critical parameter of K-SAT. For $`K=3`$, they prove $`\alpha _c4.598`$, and it is principle possible to get better bounds by including more terms in their expansion. A lower bound has also been found , $`\alpha _c>3.003`$. It is however difficult to generalize these methods to other problems, such as K-COL.

## VI The relaxation behavior

We have studied the relaxation of the energy $`ϵ`$, defined as the number of unsatisfied constraints per spin, of K-SAT and K-COL using $`T=0`$ Monte Carlo (MC) simulations and the Metropolis single-flip algorithm . A simple case where the relaxation can be understood is the ferromagnetic Ising model on a regular lattice. For this model, the energy decreases as $`ϵt^{1/3}`$ if the order parameter is conserved by the dynamics, while $`ϵt^{1/2}`$ if single spin flips that allow the total magnetization to change are used. These forms of relaxation behavior can be explained by noting that spins with the same orientation will cluster and form domains (e.g., ) with a well-defined energy. This explanation does not immediately carry over to our problems, since there is no known simple expression for the energy as a function of a length scale.

For the NPC problems studied here, the energy can not be expected always to reach zero, since there are no perfect solutions to problems with many constraints. The $`T=0`$ MC algorithm can also get stuck in local minima. A small graph where this can happen for 3-COL is shown in figure 1.

In each time step in our simulations $`N`$ spin flips are attempted. Each flip consists of selecting a random spin and randomly changing its color. The flip is allowed if this leaves the energy unchanged or lowers it, otherwise the spin is left unchanged. In temperature $`T>0`$ simulations, a flip that raises the energy $`\mathrm{\Delta }`$ units is allowed with probability $`\mathrm{exp}[\mathrm{\Delta }/T]`$.

For K-COL, we generated the problems by randomly selecting $`M`$ distinct edges from the $`\left(\begin{array}{c}N\\ 2\end{array}\right)`$ possible with $`N`$ vertices. (In graph terminology , this corresponds to using the $`𝒢(N,M)`$ model for random graphs.) For K-SAT, each clause was generated by selecting $`K`$ variables, where each variable was negated with probability $`\frac{1}{2}`$. The formula was then generated by performing this process $`M`$ times. Clauses with repeated variables were allowed in the expressions, and also repetition of clauses.

In all random systems, the question of whether or not the measured quantities are self-averaging is important. Self-averaging means that the properties of the entire (infinite) system can be understood in terms of the properties of its local subsystems. For random systems with local interactions, there is a simple argument for this (see e.g., ), but if the interactions are global, the situation is more complex. However, for spin glasses such as the Sherrington-Kirkpatrick model, the energy and other simple quantities are self-averaging, and we assume that the energy is self-averaging also in K-SAT and K-COL. This assumption of self-averaging is supported by our simulations; averaging over a small number of large graphs appears to give the same results as averaging over many small graphs. Schreiber and Martin  have recently studied various local search methods for the graph partitioning problem and found strong numerical evidence for self-averaging even for sparse graphs. They also provide some arguments for why self-averaging should hold, and conjecture that their result holds for all constraint satisfaction problems.

Our simulations show a transition between qualitatively different forms of relaxation behavior for K-SAT and K-COL. For small values of $`\gamma `$ or $`\alpha `$, we find exponential relaxation to zero energy. For larger values of the parameter, the relaxation becomes algebraic. For large enough $`\gamma `$ or $`\alpha `$, the energy freezes at a non-zero value.

The change between exponential and power law relaxation occurs for smaller values of $`\gamma `$ and $`\alpha `$ than the transition in solvability. More extensive numerical investigations will be necessary to determine the exact nature of this transition, in particular whether it is a sharp transition, or whether there are intermediate regions with different forms of relaxation behavior. In section VII below, we study the freezing transition, where the final energy of the MC algorithm changes from zero to a non-zero value. When this quantity is measured, a sharp transition is observed. This transition does not necessarily coincide with the change in functional form of the decay.

Figure 2 shows the relaxation behavior for 3-COL in a log-log diagram, while figure 3 shows that the behavior for 3-SAT is similar. The data in these figures was determined for systems of size $`10^7`$ (for 3-COL) and $`10^6`$ (for 3-SAT). The relaxation showed similar behavior for other values of $`K`$ as well.

The change between exponential and power law decay is illustrated in more detail for 3-COL in figure 4, where we have plotted $`ϵ(t)ϵ(500)`$ for systems of $`N=10^6`$ spins. The figure shows data for $`\gamma =0.5`$ up to 3.0 in increments of 0.1. For small $`\gamma `$, the decay is exponential (see also figure 5 below for $`\gamma =1`$). When $`\gamma `$ is increased, a change from exponential to power law behavior is observed. A reasonable fit to power law behavior is found approximately for $`\gamma 2`$.

Figure 5 shows the data for $`\gamma =1`$, where a reasonable fit to exponential relaxation is obtained, and for $`\gamma =1.5`$, where a cross-over behavior is seen (these and the following data were obtained for systems of size $`N=10^7`$) Figure 6 shows the data for $`\gamma =2`$, 3, and 4, where a power law $`ϵϵ_0+t^\mu `$ is found. Figure 7 shows that the power law also applies for $`\gamma =5`$ and 8; here the exponent is given by $`\mu 0.85`$.

The exponents for 3-COL are summarized in table I; note that due to finite-size effects only data up to $`t=200`$ was used to determine the exponent for $`\gamma =2`$.

The values for which the energy freezes for different $`\gamma `$ in 3-COL are shown in figure 8. For large $`\gamma `$, an approximately linear increase is observed.

For 3-SAT, $`\alpha =2`$ gives exponentially fast decay, while there appears to be a crossover behavior for $`\alpha =3`$. For $`\alpha =4`$ and 6, power law relaxation is obtained, $`ϵϵ_\mathrm{𝟎}+t^\mu `$ with $`\mu 0.6`$ in both cases, see figure 9. As in 3-COL, $`ϵ_0`$ was found to increase approximately linearly with $`\alpha `$.

No significant change in the behavior was seen in finite temperature simulations. Raising the temperature makes it possible to escape from one local minimum, but the system is then trapped in another before the ground state is reached. Raising the temperature further repeats this scenario but also increases the fluctuations in the energy. For high enough temperatures, the fluctuations take over completely and the system is not trapped in any local minimum.

We found the same form of power law relaxation with approximately identical exponents for temperatures up to $`T=0.4`$, but the frozen-in value of the energy, $`ϵ_0`$, depended on the temperature, see figure 10. For $`\gamma =4`$, the $`T=0.2`$ runs were able to achieve an almost 30% better state than the $`T=0`$ runs.

Using simulated annealing to find the ground states of some NP-complete spin glasses ($`3D`$ $`\pm J`$ and infinite range models), Grest et al  found a logarithmic decay of the energy, $`ϵϵ_0+1/\mathrm{log}t`$. Similar decay was also recently found by Kühn et al , who used an algorithm where attempts were made to flip several spins at once. The number of simultaneous flips, which plays the same rôle as the temperature in simulated annealing, was then slowly decreased. These methods will always find the ground state, whereas the $`T=0`$ MC algorithm studied here can get stuck in local minima. The faster relaxation of the $`T=0`$ MC method thus comes at the price of having no guarantee of finding the ground state.

The exact values of the critical parameter for K-COL varies depending on the ensemble of graphs used . We tried different ensembles and found no significant differences in the relaxation behavior (e.g., $`\varphi ^n`$-model Feynman diagrams show approximately the same behavior as random graphs with connectivity $`\gamma =n`$).

Since the MC algorithm can get stuck in local minima, one should verify that the behavior does not depend on the initial values of the spins. It is also important to check to what extent the result depends on the choice of random graph. To test this, we have performed simulations where in addition to averaging over several different graphs we also restarted the MC algorithm with different inital spin configurations. In particular, 3-COL with $`\gamma =2`$ and 4 was studied in detail.

In figure 11 we compare the results for $`N=10^7`$ shown above to runs for a smaller system ($`N=10^4`$), where an average over 1000 initial states was performed for each of a larger set of graphs. The figure shows the average energy from these runs and $`ϵ(t)`$ from the runs with $`N=10^7`$ for $`\gamma =2`$ and 4. A reasonable agreement is found in both cases.

The variation of the result depending on the choice of graph and initial state is further illustrated in figures 12 and 13, which show the average energy and the standard error

$`\sigma ={\displaystyle \frac{\sqrt{ϵ^2ϵ^2}}{\sqrt{N_gN_r}}},`$

(where $``$ denotes the average over graphs and restarts, and $`N_g`$ and $`N_r`$ stand for the number of different graphs and the number of runs, respectively) for $`N=10^4`$ and $`\gamma =2`$, for each of 56 different runs with $`N_r=1`$ and $`N_g=1000`$. Figure 14 shows the energy and standard error for $`N=10^4`$ and $`\gamma =4`$ for 54 different runs with $`N_r=1`$ and $`N_g=1000`$, while figure 15 shows the standard errors from each of these runs.

For $`\gamma =4`$, the variation among runs on a single graph (see figure 15) is very small compared to the variation among randomly generated graphs shown in figure 14; for $`\gamma =2`$ these are of comparable magnitude.

We also found similar relaxation behavior when we started with all spins having the same value, and did not find any differences using different random number generators. The generators used include the Mitchell-Moore additive generator (e.g., ), a 48-bit multiplicative congruential generator with multiplicative and additive constants 25214903917 and 11, respectively, and the standard C library rand() function.

## VII Comparison to other local search methods

The Monte Carlo method is an example of a local search method. Local search methods start with a candidate solution and try to improve it by changing it locally (e.g., by flipping a spin). Various other such methods have been used to study the phase transition in search cost in NPC-problems. Gent and Walsh  studied small problems using the GSAT method, a hill-climbing procedure for solving K-SAT problems, and found exponential relaxation over very small time ranges. The GSAT method is similar to MC simulations, but differs in the way a spin is chosen to flip. Instead of flipping a random spin, GSAT selects that spin whose flipping will decrease the energy the most. Clark et al  used only solvable problems (determined by first using a complete backtracking method) and found an easy-hard-easy transition in search cost for two local search methods. The hardest problems occurred approximately at the true transition.

In our case, the behavior is quite different. This is particularly evident in figure 16, which shows the relaxation for 3-COL with $`10^5`$ spins using only colorable graphs. It shows the same behavior as in figure 2 for both colorable and uncolorable problems. An exceptional case is 2-COL, where we found different behavior when only colorable graphs were used, see figure 17. The reason for this is that for $`K=2`$ and colorable graphs there is no difference between the ferromagnetic and antiferromagnetic models. Results and arguments for the relaxation behavior of ferromagnetic Potts models on random graphs will be presented elsewhere .

Problems with large connectivities are thus always hard to solve using the Monte Carlo method. This is a short-coming of the MC algorithm — even if the problem is solvable, the MC algorithm can get stuck in local minima that the other, smarter local search methods manage to avoid.

In order to quantify the difference between the MC algorithm and other local search methods, and to compare the freezing transition with the solvability transition, we determined the fraction $`\eta (t)`$ of problems for which the MC method did not find an $`ϵ=0`$ ground state in less than $`t`$ MC steps per spin. We found that $`\eta `$ displayed a behavior similar to the fraction of solvable problems — there appears to be a phase transition, but for a smaller value of $`\gamma `$.

In figure 18, we plot $`\eta (10^5)`$ for 3-COL against the rescaled parameter

$$(\frac{\gamma }{\gamma _c}1)N^{1/\nu }$$
(7)   
with $`\gamma _c=2.4`$ and $`\nu =3.75`$ and $`N`$ ranging from 20 to 1000. It is clear that there is a freezing transition well below the occurrence of the solvability transition. However, it is necessary to be cautious when drawing conclusions from small systems. The largest systems we have simulated are of size $`10^7`$ for 3-COL and $`10^6`$ for 3-SAT; these systems were used to fit the relaxation behavior shown in table I above. To determine the freezing transition, we used considerably smaller systems.

The approximate values where the transition in figure 18 was found are shown in table II for $`2K5`$, together with the experimental values for the solvability transition for K-SAT from , and values for the solvability transition for K-COL. The latter were obtained by ourselves using a non-optimized backtrack-search program with the Brelaz heuristic . For $`K=3`$, our value coincides with the literature ; we are not aware of any values for $`K=4`$ and $`5`$ in the literature. For $`K4`$, the values for the solvability transition are probably not very accurate — the results vary strongly with the number of variables and the number of graphs tested. It is clear, however, that the relaxation transition happens for smaller values than the solvability transition.

Most of our data for the freezing transition was determined by averaging over between 100 (for $`N=1000`$) and 1000 different graphs with a single MC run on each graph. We have also made some runs where we restarted the MC algorithm using different initial spin configurations for each graph. For small $`N`$, $`\eta (t)`$ varied about 13% depending on whether we used 1, 10 or 100 restarts, but for $`N80`$ there was essentially no difference between the different runs, see figure 19. In each of these runs we averaged over 1000 different graphs.

The precise location of the freezing transition depends on the algorithm used to determine $`\eta (t)`$. For example, finite temperature MC methods give a higher value for the location of the transition. This illustrates the difference between the freezing transition and the solvability transition. The former is a property of the method used to solve a problem, while the latter is a characteristic of the problem itself.

## VIII Other results

We have also found interesting scaling relations for other quantities. For graph coloring, the energy as we have defined it above measures the fraction of frustrated edges. Another possibility is to measure the fraction of frustrated spins, $`ϵ_s`$. This quantity shows the same relaxation behavior, if used in the Metropolis algorithm instead of $`ϵ`$.

It is also interesting to study the ratio $`\gamma _{\text{frust}}=\frac{2ϵ}{ϵ_s}`$, which is the connectivity of the sub-graph spanned by the frustrated edges. Figure 20 suggests a transition from $`\gamma _{\text{frust}}=1`$ for easy problems to $`\gamma _{\text{frust}}>1`$ for harder. If the connectivity is 1, all frustrated edges are isolated, while connectivity 2 would mean that all frustrated spins were connected in chains.

The Monte Carlo dynamics itself has interesting properties. Figure 21 shows the fraction $`r(t)`$ of persistent spins, i.e., those that have not yet been flipped, as a function of time for 3-SAT. The data suggest a transition from an exponentially fast to a logarithmically slow decay as $`\alpha `$ is increased (this will be explored further elsewhere). For the Ising and Potts models on a square lattice, this quantity has been found to scale with a power law dependence on time .

## IX Heuristic arguments for power law relaxation

In K-COL and K-SAT there are global barriers to local improvement. In order to lower the frustration for a spin, we have to make changes to many other, unfrustrated spins (compare figure 1). One way of explaining the relaxation is to look at these barriers and see how long it takes to overcome then.

In the following, we call a proposed assignment of the $`N`$ variables a state in the space of all solutions. Consider the probability $`p(i)`$ that a state has $`i`$ unsatisfied constraints, so that $`ϵ=i/N`$. In analogy with expression (1) for K-COL, this can be approximated by

$$p_{\text{sol}}(i)=(1p_v)^{Mi}p_v^i\left(\begin{array}{c}M\\ i\end{array}\right),$$
(8)   
where $`p_v`$ is the probability that a constraint is violated. The binomial factor $`\left(\begin{array}{c}M\\ i\end{array}\right)`$ represents the number of different ways to choose the $`i`$ unsatisfied edges. Recall that for K-COL $`p_v=1/K`$ and for K-SAT $`p_v=1/K^2`$. In the approximation, we neglect all correlations between the constraints, such as the presence of triangles and other regular structures in the graph to be colored. In the following, only K-COL will be considered.

The Monte Carlo algorithm works by generating a new state, and changing to this state if its energy is lower than or equal to that of the current state. In the approximation we assume that the probability that the new state has energy $`e^{}`$ is proportional to $`p_{\text{sol}}(e^{})`$. Given a state with energy $`e`$, all transitions to states with energy $`e^{}>e`$ are forbidden. The probability of staying in a state with energy $`e`$ is then proportional to $`_{i=e}^Mp_{\text{sol}}(i)`$, while the probability for a transition to an energy $`e^{\prime \prime }<e`$ is proportional to $`p_{\text{sol}}(e^{\prime \prime })`$.

This process can be viewed as a Markov chain. The transition matrix between states of different energy has elements $`p_{ij}`$ given by the probability that state $`i`$ is followed by state $`j`$:

$$p_{ij}=\{\begin{array}{cc}0,& \hfill i<j\\ _{n=i}^Mp_{\text{sol}}(n)=1_{n=0}^{i1}p_{\text{sol}}(n),& \hfill i=j\\ p_{\text{sol}}(j),& \hfill i>j\end{array}$$
(9)   
for $`i,j=0\mathrm{}M`$. The transition probability $`p_{ii}`$ can be written in terms of the hypergeometric function $`{}_{2}{}^{}_{1}^{}(a,b,c;z)`$  as

$`p_{ii}=p_{\text{sol}}(i)_2_1(1,iM,i+1;{\displaystyle \frac{p_v}{p_v1}}).`$

This approximation ignores the detailed mechanisms of the MC-algorithm (which can only change a state locally), instead we consider transitions between classes of states with the same energy.

Using this transition matrix to evolve the states in time numerically, power law relaxation behavior is found, see figure 22. The exponent differs from that in our MC simulations. One reason for this is that the MC simulations only allow single-spin flips, while equation (9) allows jumps between arbitrary spin configurations — the global and the local energy landscapes are different. We have performed simulations using a MC algorithm that changes the entire spin configuration instead of just a single spin. The results indicate a slower power law relaxation, with an exponent that is consistent with figure 22. In contrast to the regular MC simulations, the exponent here showed a weak dependence on system size. The dependence was the same for both the simulations and time-evolution using equation (9).

In the MC algorithm, $`N`$ attempts to flip a spin are made in each time step. Inspired by statistics of the energy landscape, we have made simulations where we require that there are $`N`$ accepted spin flips in each time step (in $`T=0`$ simulations, this means that time is not increased for flip attempts to higher energy.) Time was increased by $`1/N`$ for each change in the spin configuration. The power law relaxation behavior did not change when this modified algorithm was used.

For these simulations, we can ignore all transitions to states with higher energy. Since a single spin flip is unlikely to decrease the energy by more than one (see figure 20), we can also ignore transitions to states with energy $`j<i1`$, and consider a two-state system.

To estimate the time needed to go from state $`i`$ to $`j=i1`$, we assume that $`p(i)/p(i1)`$ can be approximated by

$$p_{\text{sol}}(i)/p_{\text{sol}}(i1)=\frac{p_v}{1p_v}\frac{Mi+1}{i}\frac{1}{ϵ}.$$
(10)   
Since the average time needed to pass from state $`i`$ to state $`j=i1`$ is given by the inverse of the transition probability per unit time, we then have

$$\mathrm{\Delta }t=\frac{p(i)+p(i1)}{p(i1)}=1+\frac{p_{\text{sol}}(i)}{p_{\text{sol}}(i1)}1+\frac{1}{ϵ}$$
(11)   
for the time needed to go from state $`i`$ to $`i1`$.

Another viewpoint is to consider the free energy barriers, $`\mathrm{\Delta }f`$, in the problem. Since transitions to states with higher energy are not allowed by the $`T=0`$ Monte Carlo dynamics, we can neglect the energy part of $`\mathrm{\Delta }f`$. The entropic part will however be positive,

$$\mathrm{\Delta }f=k_bT\mathrm{log}N(j)+k_bT\mathrm{log}N(i)=k_bT\mathrm{log}\frac{p_{\text{sol}}(i)}{p_{\text{sol}}(j)}$$
(12)   
for a transition from $`i`$ to $`j`$. The time needed to overcome this free energy barrier is

$$\tau \mathrm{exp}(\frac{\mathrm{\Delta }f}{k_bT})\frac{p_{\text{sol}}(i)}{p_{\text{sol}}(j)}\frac{1}{ϵ},$$
(13)   
if $`i=j+1`$ and $`ϵ=i/N`$.

By inverting equations (11) and (13), power law relaxation with an asymptotic $`t^1`$ dependence can be obtained. This is justified for large $`\gamma `$, where all free energy barriers have roughly the same size, but not for $`\gamma `$ below the freezing transition. Similar difficulties also arise, e.g., for coarsening, where arguments for relaxation usually are given for nucleation, where a single domain grows or shrinks. It is then justifiable to invert the equation relating the energy of the domain and the time needed to shrink it. In the coalescence regime, however, several such domains form and grow separately before they coalesce. The inversion is then not allowed, even though the same relaxation behavior is found in simulations.

The power law with $`t^1`$ dependence approximately matches the behavior of K-COL for critically constrained $`\gamma `$’s (see table I). For overconstrained problems, the relaxation is slower than $`t^1`$. This could depend on the local energy landscape being different from the global when the problems are overconstrained. The arguments above should be viewed only as a simplistic first approximation, and cannot be expected to reproduce all important features of the numerical results.

The barriers in equation (12), which are proportional to the logarithm of the energy, do not appear in any of the coarsening classes introduced by Lai et al . The difference seems to be that here there are entropy barriers. That the entropy is high for low energies is explained by the fact that the interactions are antiferromagnetic and the variables are Potts spins. If each spin can have $`K`$ different values and the interactions are ferromagnetic, each satisfied bond can be satisfied in $`K`$ different ways, while there are $`K^2K`$ choices for the unsatisfied ones. If the interactions are antiferromagnetic, each satisfied bond has $`K^2K`$ possibilities and the unsatisfied $`K`$. Since there are more satisfied bonds than unsatisfied, an antiferromagnetic problem is less restricted, at least for $`K>2`$. As always, this reasoning ignores all correlations between constraints in the problem, such as triangles in graphs.

## X Conclusions and Discussion

We have shown that there are additional transitions in the behavior of two NPC problems — graph colorability and satisfiability. There is a freezing transition, and the relaxation changes from exponential to power law as the difficulty of the problem increases. These transitions occur for smaller values of the parameters $`\gamma `$ and $`\alpha `$ than the transitions in solvability and search cost.

We also studied the effects of using only colorable graphs, and measured other quantities, such as the fraction of persistent spins for K-SAT, which also showed a transition.

A simple heuristic argument for the form of the relaxation behavior depending on entropy barriers was also given. The free energy barriers in these problems have a less pronounced energy dependence than in other models. One can compare them to systems where logarithmic relaxation has been found for $`T0`$ and $`T=0`$ MC simulations, such as in , where a tiling model containing no randomness at all was studied. The main difference between these systems and ours seems to be that they have local interactions and frustration, whereas we have infinite-range interactions and frustration.

The problems with logarithmic relaxation are in some ways simpler than the NPC problems studied here. Their ground states can always be found in polynomial time, and they do not contain any difficult optimization problem. On the other hand, K-SAT and K-COL have ground states that may require exponential time to find. Our results show that it can be easier to find approximate solutions to these hard problems than for the easy problems giving logarithmic relaxation. Of course, to actually find the ground state with the MC algorithm is very difficult for K-COL and K-SAT. What the NPC problems gain in finding an approximate solution is lost when it comes to finding the ideal, ground state solution.

Other problems where entropy barriers and logarithmic relaxation are found include the backgammon model introduced by Ritort  and the random-field Ising model . The backgammon model contains no energy barriers at all, only entropy barriers, and can be solved  by considering random-walk models with entropy barriers. In our models, the entropy barriers scale as the logarithm of the energy. This means that the time needed to overcome the free energy barrier loses its exponential dependence and becomes algebraic, leading to power law relaxation. We have also found this behavior in more general constraint satisfaction problems than K-COL and K-SAT . Relaxation behavior similar to that in our models has also been obtained by Campellone et al for a short-range $`p`$-spin glass model . These models should be studied further.

We see several possibilities for future work in this field. More detailed studies of the behavior of the fraction of persistent spins as well as of damage spreading and hysteresis in strongly disordered systems are needed. Since hard optimization problems are also spin glasses, a search for spin glass characteristics such as aging and chaotic responses to small perturbations in the hamiltonian may also be interesting. We are preparing a study of the energy landscapes of these and related models. Investigating the connection between the topology of the random graph and that of the energy landscape may also provide more insight into why NPC problems are difficult.

## Acknowledgements

We would like to thank an anonymous referee for useful comments.