Sudoku as k-SAT

Because our continuous-time dynamical system15 was designed to solve k-SAT formulae in conjunctive normal form (CNF), we first briefly describe how Sudoku can be interpreted as a +1-in-9-SAT formula and then how it is transformed into the standard CNF form. Further details are shown in the Methods section.

In a Sudoku puzzle we are given a square grid with 9 × 9 = 81 cells, each to be filled with one of nine symbols (digits) , i, j = 1, …, 9 (with the upper-left corner of the puzzle corresponding to i = 1, j = 1). When the puzzle is completed, each of the columns, rows and 3 × 3 sub-grids (blocks partitioned by bold lines, Fig. 1a) must contain all the 9 symbols. Equivalently, all 9 symbols must appear once and only once in each row, column and 3 × 3 sub-grid.

Figure 1 Sudoku and its boolean representation. (a) a typical puzzle with bold digits as clues (givens). (b) Setup of the boolean representation in a 9 × 9 × 9 grid. (c) Layer L 4 of the puzzle (the one containing the digit 4) with 1-s in the location of the clues and the regions blocked out for digit 4 by the presence of the clues (shaded area). Full size image

To formulate Sudoku as a constraint satisfaction problem (CSP) using boolean variables, we associate to each symbol (digit) an ordered set of 9 boolean variables (TRUE = “1”, FALSE = “0”). The digit D ij in cell (i, j) will be represented as the ordered set ( ) with , a = 1, …, 9, such that always one and only one of them is 1 (TRUE). In particular, D ij = a is set to be equivalent to , where δ a,b is the Kronecker delta function. This way we have in total 9 × 9 × 9 = 729 boolean variables , which we can picture as being placed on a 3D grid (Fig. 1b), with a corresponding to the grid index along the vertical direction and hence a is the digit that is filling the corresponding (i, j) cell in the original puzzle. The corresponding horizontal 9 × 9 2D layer at height a will be denoted by L a . Introducing the notion of such (horizontal) layers makes it easier to express the constraints of the Sudoku rules on its representation by 0-s and 1-s as described below. For example in the puzzle shown in Fig. 1a, we have D 1,9 = 4. In the given vertical column the variable in the ath cell (that is in layer L a ) is (shown as the boolean variable 1 filling the cell next to the digit 4 shown in red, in Fig 1b). This setup allows us to encode the Sudoku constraints in a simple manner. They come from: 1) uniqueness of the digits in all the (i, j) Sudoku cells, 2) a digit must occur once and only once in each row, column and in each of the nine 3 × 3 subgrids and 3) obeying the clues. Constraint type 1) was already expressed above, namely that for every cell (i, j), in the set ( ) one and only one variable is TRUE, all others must be FALSE. Type 2) constraints are similar, e.g., in row i and layer L a the set ( ) must contain one and only one TRUE variable, all others must be false and this must hold for all rows and layers, etc. Observe that all constraints are in the form of a set of 9 boolean variables of which we demand that one and only one of them be TRUE, all others FALSE. When this is satisfied, we say that the constraint itself (or “clause”) is satisfied, or TRUE. Such CSPs are called +1-in-k-SAT and they are part of so-called “locked occupation problems”, which is a class of exceptionally hard CSPs18,19. Type 3) constraints are generated by the clues (or givens) which are symbols already filled in some of the cells and their number and positioning determines the difficulty of the puzzle. They are also set in a way to guarantee a unique solution to the whole puzzle. If there are given d clues, then this implies setting d boolean variables to TRUE, which means eliminating exactly 4d constraints of type 1) and 2) (one vertical or uniqueness constraint, one row, one column and one 3 × 3 subgrid constraint). Thus, Sudoku is a +1-in-9-SAT type CSP with N boolean variables and 324 − 4d constraints of +1-in-k-SAT type (k ≤ 9). N is a complicated function of the positioning of the clues. The example in Fig. 1a has d = 22 clues with N = 232 unknown boolean variables. In layer L 4 as illustrated in Fig. 1c there are 28 unknown boolean variables (white cells). These 28 variables appear in a total of 17 constraints of +1-in-k-SAT type. More precisely there is one +1-in-2-SAT, six +1-in-4-SAT, four +1-in-5-SAT and six +1-in-6-SAT type of constraints related to L 4 . The other layers generate the remaining 324 − 4d − 17 = 219 of +1-in-k-SAT type constraints (with k ≤ 9).

Since our continuous-time SAT solver has been designed to solve boolean satisfiability problems in conjunctive normal form (CNF), we need to bring the +1-in-k-SAT type problems above into this form and thus formulate it as a k-SAT CNF problem. The CNF is a conjunction (AND, denoted by ∧) of clauses each clause expressed as the disjunction (OR, denoted by ∨) of literals. For k-SAT in CNF there are N boolean variables x i = {0, 1} and an instance is given as a propositional formula , which is the conjunction of M clauses C m : , with , k m ≤ k and z j is a literal, that is ( is the negation of x j ). According to a well known theorem of propositional calculus, all boolean propositions can be formulated in CNF using De Morgan's laws and the distributive law and thus any +1-in-k-SAT type constraint as well. The Methods section describes how to translate +1-in-k-SAT type constraints into CNF.

Once the transformation to CNF is completed we are left with N variables and M clauses of the type described above, called CNF clauses from here on. In our case the number of variables appearing in a CNF constraint has the property 1 ≤ k m ≤ 9. The parameters N, M and depend on the clues that are difficult to express analytically, but easy to determine computationally. The puzzle from Fig. 1a is ultimately formulated as a CNF SAT problem with N = 232 variables and M = 1718 CNF clauses. An often used parameter of a satisfiability problem is the number of CNF constraints per variable, or constraint density, α = M/N, also used as a typical hardness indicator, however, as we show below, this is not an accurate measure of hardness.

The continuous-time deterministic k-SAT solver

In Ref [15] a continuous-time deterministic solver was introduced to solve k-SAT problems in conjunctive normal form. The set of clauses specifying the constraints are translated into an M × N matrix: C = {c mi } with c mi = 1 if the variable x i is present in clause m in direct (non-negated) form, namely , c mi = −1 if and c mi = 0 if x i and are both absent from C m . To every variable x i one associates a continuous spin variable such that when s i = ±1 then and to every clause C m one associates the function:

We have for all . It is easy to check that K m = 0 only for those values for which the corresponding x i -s satisfy clause C m (otherwise we always have K m > 0). That is, K m plays the role of an energy function for clause C m and its ground state value of K m = 0 is reached if C m is TRUE and only then. We also need the quantities that is, with the i-th term missing from the product in (1). Clearly, . The continuous time dynamical system introduced in15 is defined via the set of (N + M) ordinary differential equations (ODEs):

with the only requirements that , ∀i and a m (0) > 0, ∀m. The latter implies from (3) that a m (t) > 0, ∀m, t. It was shown in Ref [15] that system (2–3) always finds the solutions to k-SAT problems (encoded via the C matrix), when they exist, from almost all initial conditions (the exception being a set of Lebesgue measure zero). Here we give an intuitive picture for why that is the case. Due to (3) the auxiliary variables a m grow exponentially at rate K m . That is, the further is K m from its ground-state value of 0, the faster a m grows (in that instant). Moreover, the longer has K m been away from zero, the larger is a m , as seen from the formal solution to (3): . Equation (2) can equivalently be written as a gradient descent on an energy landscape V(s, a), that is ds/dt = − ∇ s V, where ∇ s is the gradient operator in the spin variables and . Clearly, V ≥ 0 ∀t and V = 0 if and only if s is a k-SAT solution, i.e., satisfies all the clauses (K m (s) = 0, ∀m).

From the behavior of the a m variables discussed above it also follows that the least satisfied constraints will dominate V (terms with the largest a m -s). Without restricting generality, let the term be the most dominant at t. Then keeping only the dominant term on the rhs of (2) for those i for which c 1i ≠ 0 we get or, equivalently: . This shows that the term (1 − c 1i s i ) is driven exponentially fast towards zero, that is towards satisfying K 1 (and all the other constraints containing this term). As K 1 decreases, some other constraint becomes dominant and thus, in a continuous fashion, all constraints are driven towards satisfiability. The exponential growth guarantees that the trajectory is always pulled out of any potential well. When the problem is unsatisfiable, the system generates a chaotic dynamics in [−1, 1]N, indefinitely. For more details about the properties of the continuous-time dynamical system (CTDS) (2–3) see Ref. [15].

Puzzle hardness as transient chaotic dynamics

Since Sudoku puzzles always have a solution, the corresponding boolean SAT CNF formulation also has a solution and system (2-3) will always find it. The nature of the dynamics, however, will depend on the hardness of the puzzle as we describe next.

In Fig.2a we show an easy puzzle with 34 clues (black numbers)20. After transforming this problem into SAT CNF, we obtain N = 126 and M = 717, with a constraint density of α = M/N = 5.69. As described above, in our implementation there is a continuous spin variable associated to every boolean variable in every 3D cell (i, j, a). In the right panels of Fig. 2 we show the dynamics of the spin variables in the cells of the 3 × 3 grid formed by rows 4–6 and columns 7–9. The curves are colored by the digit a they represent (a = 1, …, 9) as indicated in the color legend of Fig 2. The dynamics was started from a random initial condition. Indeed, our solver finds the solution very quickly, in about 15 time units, for the easy puzzle in Fig.2a.

Figure 2 Solving Sudoku puzzles with the deterministic continuous-time solver (2–3). (a) presents an easy puzzle with the evolution of the continuous-time dynamics shown within a 3 × 3 grid (rows 4–6, columns 7–9). (b) shows the same, but for a known, ultra-hard puzzle called “Platinum Blonde”21,22,23. The trajectories in the right panels show the evolution of the analog variables colored by the corresponding digit a. Thus for each cell (i, j) we have 9 such running trajectories, but they cannot always be discerned as many of them are running on top of each other, close to -1, as seen in (a). Full size image

In Fig. 2b we show the dynamical evolution of variables for an ultra-hard Sudoku instance with only 21 clues. This puzzle has been listed as one of the world's hardest Sudokus and even has a special name: “Platinum Blonde”21,22,23 and it was the most “difficult” for our solver among all the puzzles we tried. After transforming it into SAT CNF, we obtain N = 257 variables and M = 2085 constraints. Not only that we have twice as many unknown variables but the constraint density α = M/N = 8.11 is also larger than in the previous case, signaling the hardness of the corresponding SAT instance. The complexity of the dynamics in this case is seen in the right panel of Fig. 2b, exhibiting long chaotic transients before the solution is found at around . For an animation of the dynamics for a similarly hard puzzle14 see Ref [24].

We can also observe from the right panels in Fig 2 that there is one dominating digit (a-value), corresponding to which vertical cell at that given (i, j) grid cell has the largest value. This can be taken as the digit D ij the solver is considering in the given grid cell (i, j) at that moment. We will use this observation to provide below an alternate illustration of the dynamics' transiently chaotic behavior. Let us fix a random initial condition except for two chosen variables that are varied along the points of a square grid within the domain [−1, 1]2. There is no particular relevance as to which pairs of variables are chosen to be varied, let us denote them by s 1 and s 2 . Let us choose an arbitrary empty cell (p, q) in the original Sudoku puzzle and monitor the dominating digit in it at time t. We will color the initial conditions in the plane (s 1 , s 2 ) according to the dominating digit in (p, q) at time t. This will provide a map expressing the “sensitivity to initial conditions” that varies across time. Since all puzzles have solutions, the maps eventually assume one solid color according to the digit of the solution in the monitored cell, however, for hard puzzles, it may assume highly complex patterns before it does that, as shown in Fig. 3. Here we show these colormaps for the easy and hard Sudoku puzzles from Fig. 2 at times t = 10, 15, 20. For the easy puzzle (top row) the cell was chosen to be (p, q) = (1, 1). At time t = 10 the whole map shows D 1,1 = 6 (orange), which is not the solution digit (it is still searching for the solution). At time t = 15, however, we see two clearly separated domains, in one of them D 11 = 6, in the other D 11 = 9 (light blue), the latter being the correct digit. As time passes, the orange (incorrect) domain shrinks, because trajectories from an increasing number of initial conditions find the solution. At t = 20 almost the whole map shows the correct digit D 11 = 9, except for a thin line.

Figure 3 Puzzle hardness as chaotic dynamics. We color the points of a 103 × 103 grid in an arbitrary plane (s 1 , s 2 ) at time instant t according to the digit D pq the solver is considering in an arbitrary but fixed cell (p, q) at that instant, given that we started the trajectory of the CTDS from those grid-points. For these initial conditions only the points in the (s 1 , s 2 ) plane were varied, all other spin values were kept fixed at the same randomly chosen values. For the same easy problem as in Fig. 2 and for (p, q) = (1, 1) (top row of panels) almost all initial conditions in this plane involve only two digits and after t = 20 the corresponding trajectories have converged to the solution digit (9, light blue), except for a thin line, which, however, will also become light blue. The bottom row of panels shows the same for the hard problem of Fig. 2 based on what happens in the cell (p, q) = (6, 8). The strong sensitivity to initial conditions appears as fractal structures of increasing complexity as time goes on, before eventually everything converges to the same color/digit (dark blue, corresponding to digit 4, not shown). Full size image

In the case of the hard Sudoku puzzle (bottom row in Fig. 3, (p, q) = (6, 8)) more colors enter the picture over time, in a fractal-like pattern. Eventually the whole map becomes dark blue (digit 4) corresponding to the solution in this cell. Changing the initial condition slightly about this fractal set (which is really a time-dependent fractal boundary) may result in a completely different digit (color) being considered in cell (p, q) at time t. This sensitivity to initial conditions is indicative of the chaotic behavior of the (deterministic) search dynamics.

The appearance of transient chaos is a fundamental feature of the search dynamics and can be used to separate problems by their hardness. In Ref [15] we have shown that within the thermodynamic limit (N → ∞, M → ∞, α = M/N = const.) of random k-SAT ensembles this appears as a phase transition at the so-called chaotic transition point α χ in terms of the constraint density α = M/N. Since there is no “thermodynamic limit” for 9 × 9 Sudoku problems (N < 729), one cannot define a simple order-parameter and use it to rate problem hardness in the same way15. However, once a problem is given, the corresponding dynamical system (2–3) is well defined and so is its dynamical behavior. Even though we do not have a well-defined ensemble-based statistical order parameter (which has little meaning for specific SAT instances anyway), here we show next how can we use a well-known invariant quantity from non-linear dynamical system's theory to categorize problem hardness for specific instances.

A Richter-type scale for Sudoku hardness

As suggested by the two examples in Fig. 2, the hardness of Sudoku puzzles correlates with the length of chaotic transients. A consistent way to characterize these chaotic transients is to plot the distribution of their lifetime. Starting trajectories from many random initial conditions, let p(t) indicate the probability that the dynamics has not found the solution by time t. A characteristic property of transient chaos11,12 in hyperbolic dynamical systems is that p(t) shows an asymptotic exponential decay: p(t) ~ e–κt, where κ is called the escape rate. For an easy to read text on transient chaos see the book13. The escape rate, a measurable quantity, theoretically can be expressed as a zero of the spectral determinant of the evolution operator corresponding to the dynamical system (2–3) and well approximated using the machinery of cycle expansions based on dynamical zeta functions12. It is an invariant measure of the dynamics in the sense that it characterizes solely the chaotic nonattracting set in the phase space of the system and it does not depend on the distribution of the initial conditions, its support, or the details of the region from where the escape is measured (as long as it contains the non-attracting set)11.

In Fig. 4a we plot the distribution p(t) in log-linear scale for several puzzles gathered from the literature. The distributions were obtained from over 104 random initial conditions. The decay shows a wide range of variation between the puzzles. For easy puzzles the transients are very short, p(t) decays fast resulting in large escape rates but for hard puzzles κ can be very small. Fig. 4b shows a zoom onto the p(t) of hard puzzles. For example, for the puzzle in Fig.2a we obtain κ = 0.156, whereas for Fig.2b (Platinum Blonde) the escape rate is κ = 0.00026. In spite of the large variability of the decay rates, we see that in all cases the escape is exponentially fast (the faster than exponential appearance for some of the curves in Figures 4a,b is due to finite size statistics and the finiteness of the time interval considered).

Figure 4 Escape rate as hardness indicator. (a) shows the distribution in log-linear scale of the fraction p(t) of 104 randomly started trajectories of (2–3) that have not yet found a solution by analog time t for a number of Sudoku puzzles taken from the literature (see legend and text) with a wide range of human difficulty ratings. The escape rate is obtained from the best fit to the tail of the distributions. (b) is a magnification of (a) for hard puzzles. (c) and (d) show the escape rate κ in semilog scale vs the number of clues d and constraint density α indicating good correlations with human ratings (color bands). (e) shows the relationship between the number of clues d and α for the puzzles considered. Full size image

The several orders of magnitude variability of κ naturally suggests the use of a logarithmic measure of κ for puzzle hardness, see Fig.4c, which shows the escape rates on a semilog scale as function of the number of clues, d. Thus, the escape rate can be used to define a kind of “Richter”-type scale for Sudoku hardness:

with easy puzzles falling in the range 0 < η ≤ 1, medium ones in 1 < η ≤ 2, hard ones in 2 < η ≤ 3 and for ultra-hard puzzles η > 3.

We chose several instances from the “Sudoku of the Day” website20 in four of the categories defined there: easy (black square), medium (red circle), hard (green ×) and absurd (blue star). These ratings on the website try to estimate the hardness of puzzles when solved by humans. These ratings correlate very well with our hardness measure η, giving an average hardness value of 〈η〉 = 0.816 for easy, 〈η〉 = 1.439 medium, 〈η〉 = 1.782 for hard and 〈η〉 = 1.809 for what they call absurd. Another site we analyzed puzzles from is “Extreme Sudoku”25 (brown + signs on Fig.4). It claims to offer extremely hard Sudoku puzzles, their categories being: evil, excessive, egregious, excruciating and extreme. Indeed those puzzles are difficult with a range of η ∈ [1.1, 1.9] on the hardness scale, however, still far from the hardest puzzles we have found in the literature. Occasionally, daily newspapers present puzzles claimed to be the hardest Sudoku puzzles of the year. In particular, the escape rate for the Caveman Circus 2009 winner26 (turquoise diamond) and the Guardian 2010 hardest puzzle27 (maroon diamond) are indeed one order of magnitude smaller than the hardest puzzles on the daily Sudoku websites, placing them at η = 2.93 and η = 2.82 on the hardness scale. The USA Today 2006 hardest puzzle28, however, does not seem to be that hard for our algorithm having η = 2.17 (magenta diamond). Eppstein29 gives two Sudoku examples (orange left-pointing triangles) while describing his algorithm, one with η = 1.288 and a much harder one with η = 2.017. Elser et al.14 present an extremely hard Sudoku (black filled circle), which has an escape rate of κ = 0.0023 resulting in η = 2.639.

The smallest escape rates we have found are for the Sudokus listed as the hardest on Wikipedia22,23 (red triangles). The five puzzles, which we tested are called Platinum Blonde, Golden Nugget, Red Dwarf, coly013 and tarx0134. They have a hardness in the range 3 < η < 3.6, the ultra-hard Platinum Blonde (shown in Fig.2b) being the hardest with η = 3.5789 (κ = 0.00026).

While the escape rate correlates surprisingly well with human ratings of Sudoku hardness, it is natural to expect a correlation with the number of clues, d. Indeed, as a general rule of thumb, the fewer clues are given, the harder the puzzle, however, this is not universally true1. Here we tested a few instances with minimal30, that is 17 clues and almost minimal 18 clues (orange filled circles)31,32,33. As seen from Fig.4c, these are actually easier (1.2 < η < 2.4) than the hardest instances with more d = 21, 22 clues. In Fig.4d we then plot the escape rate as function of the constraint density α = M/N, leading to practically the same conclusion. This is because the constraint density α is essentially linearly correlated with the number of givens d, as shown in Fig.4e. The apparent non-monotonic behavior of puzzle hardness with the number of givens or constraint density is due to the fact that hardness cannot simply be characterized by a global, static variable such as d or α, but it also depends on the positioning pattern of the clues1.