{'fraction_non_alphanumeric': 0.026183172655565294, 'fraction_numerical': 0.008709465381244522, 'mean_word_length': 4.5763591936469155, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Difficult Sudoku Puzzles Created by Replica Exchange Monte Carlo Method

An algorithm to create difficult Sudoku puzzles is proposed. An Ising spin-glass like Hamiltonian describing difficulty of puzzles is defined, and difficult puzzles are created by minimizing the energy of the Hamiltonian. We adopt the replica exchange Monte Carlo method with simultaneous temperature adjustments to search lower energy states efficiently, and we succeed in creating a puzzle which is the world hardest ever created in our definition, to our best knowledge. (Added on Mar. 11, the created puzzle can be solved easily by hand. Our definition of the difficulty is inappropriate.)


I. INTRODUCTION
Sudoku, which is also called Number Place, is a kind of pencil puzzles [1]. Each Sudoku puzzle has 9 × 9 grid, which has nine 3 × 3 subgrids. There are some numbers in cells. The objective of the puzzle is to complete the grid by filling numbers in empty cells so that each row, column, and subgrid contains all of numbers from 1 to 9. Sudoku puzzles are now extremely popular in the world, and recently Sudoku have attracted much attention as mathematical and physical points of view. In 2002, solving Sudoku puzzle is proved to be NP-complete problems [2]. All possible Sudoku solutions is enumerated by Felgenhauer and Jarvis [3]. Solving Sudoku puzzles corresponds to find the grand state of the antiferromagnetic 9-state Potts model with special interactions. The similarity between the Sudoku problems and spin-glass systems has been pointed out [4]. Williams and Ackland defined a Sudoku Hamiltonian, and observed thermodynamic phase transitions by utilizing Monte Carlo (MC) simulations [5]. They also pointed out that the energy landscape of the Sudoku Hamiltonian is rugged, and the model show similar behavior to spin-glass systems. From the view point of the computing, it is rather easy to find the solution of the given puzzle. However, it is not trivial to make difficult puzzles by using computers. In 2010, Dr. Arto Inkara created a difficult Sudoku puzzle (Inkara2010) [6] which is shown in Fig. 1 (a). Later, he created a more difficult one in 2012 (Inkara2012) [7] which is shown in Fig. 1 (b). To our best knowledge, Inkara2012 is the world hardest Sudoku puzzle ever created. The purpose of the present manuscript is to create a Sudoku puzzle which is more difficult than Inkara2012, and consequently is the world hardest, by utilizing a MC method.
This manuscript is organized as follows. A method to solve Sudoku puzzles and a definition of difficulty are described in Sec. II. The algorithms for finding difficult puzzles are explained in Sec. III. The numerical results are given in Sec. IV and a summary and a discussion of further issues are given in Sec. V.

II. DEFINITION OF DIFFICULTY
Difficulty of Sudoku puzzles depends on a method to solve them. Therefore, we describe the method to solve Sudoku puzzles. While there are many kinds of techniques, we adopt only two of them, pencil marks and recursive backtracking (see Fig. 2).
Pencil Marks: Pick up an empty cell. Check all the numbers in the row, column, and subblock to which the cell belongs. Then list up all numbers which are still possible in the empty cell. These numbers are called pencil marks. If a cell has only one pencil mark, then the mark is the value of the cell. Repeat the above procedure until all empty cells have two or more pencil marks.
Recursive Backtracking: Pick up the cell which has the smallest number of pencil marks. Then choose one of the pencil marks and assume that it is the value of the cell, and continue to solve the problem recursively. If the assumed value does not lead to the complete solution, then choose another value of the pencil marks.
The algorithm to solve Sudoku puzzles involves the iteration of pencil marks and recursive backtracking. Since the computational costs of pencil marks are cheap, we do not includes pencil marks to difficulty. We define three kinds of properties for Sudoku puzzles, depth, normal width, and average width, respectively. The depth is the smallest number of recursive backtracking to solve a puzzle. If a puzzle is solved only with pencil marks, then the depth of the puzzle is zero. The width is the number of recursive backtracking required to confirm that the puzzle has a unique solution. It is possible that there are two or more empty cells which have the same number of pencil marks. In Fig. 2, the smallest number of pencil marks in cells is two, and there are four cells which have two pencil marks. A value of width depends on a choice of cells to perform recursive backtracking. We define two kinds of width, the normal width and the average width. The recursive backtracking process constructs a tree structure. The root node of the tree is the given puzzle. When a empty cell is chosen for recursive backtracking, then the number of pencil marks in the cell corresponds to the number of edges of the nodes. Each child node corresponds to a grid with an assumed number in the chosen cell. The tree is constructed recursively for children. Since we have a choice of cells to perform recursive backtracking, there are many possibility for the tree structure. The depth is the shortest path from the root node to the node describing answer between all possible trees. The width is a number of nodes in the tree graph. Computational cost to search the solution is proportional to width of the graph. Therefore, we adopt width as difficulty of Sudoku puzzles. Since the value of the normal width depends on orientation of a puzzle, the average width is more appropriate for difficulty of Sudoku puzzles than the normal width. Since a number of possible trees increases exponentially as a number of recursive backtracking increases, it is impossible to enumerate all possible trees. Therefore, we estimate the average width by MC sampling. In the present manuscript, we calculate average width for each puzzle from 100 independent samples.

A. Model
Suppose a solution of a puzzle, i.e., all of cells are filled with numbers, is given. To create a puzzle, we have to remove some of numbers from the grid. The rest numbers are hints of this puzzle. We label cells from 1 to 81, and describe the state of i-th cell by a spin s i ; the i-th cell is empty when s i = 0 and the cell keeps the number of the answer when s i = 1. A set of spin configuration {s i } denotes a puzzle. We define an Ising spin-glass like Hamiltonian of this system as where U denotes the internal energy given by the configuration {s i }, J(> 0) is the interaction energy, and h(> 0) is the amplitude of the external field, respectively. The interaction energy U is defined so that the energy of the system decreases when the puzzle is more difficult. We define two kinds of internal energy, depth energy U d and width energy U w . Consider a puzzle which depth is d and normal width is w. Then the depth and width energies are defined to be Since J > 0, the energies decrease as depth or width increase. We adopt logarithm for the width energy since width increases exponentially as a number of recursive backtracking increases. The energies depend on spin configuration {s i }, but the relation between them is highly complicated. While one can easily calculate the energies from the given spin configuration, it is almost impossible to find the grand state, i.e., to find the spin configuration which describes the most difficult puzzle for the given solution of the puzzle [8]. This property is similar to that of spin-glass models. The second term in the right-hand side of Eq. (1) increase energy when a number of hints of a puzzle increases. This term corresponds to an external field.
While we have to minimize only the internal energy for the purpose to search a difficult puzzle, we added this term since to apply a bias to reduce a number of hints. Since the spins are all-up in the initial states, it is difficult to decrease energy without this term. Our purpose is to create a puzzle width a large value of average width. But we find that it is difficult to obtain large values of width from the initial state where all spins are up. Therefore, we first perform MC simulations in the depth-first order, then we switch to that in the width-first order.
In order to update a configuration of spins, we adopt the Markov Chain Monte Carlo method. Choose a spin randomly and flip it with the Metropolis criterion with the Boltzmann weight [9]. Note that, a spin-flip from down to up always increases energy, and that from up to down always decreases energy. If a new configuration has two or more solutions, then the trial is rejected.

B. Adjusting Temperature Set
Since the model is similar to Ising spin-glass models, there are many local minima in the energy land scape. In order to search lower-energy states efficiently, we adopt replica-exchange Monte Carlo (REMC) method [10,11], which is also called Parallel Tempering method. The REMC method is proposed by Hukushima and Nemoto in order to study hardly-relaxing systems such as spinglass systems efficiently. In the REMC method, many replicas sharing the identical Hamiltonian are simulated simultaneously and temperatures of replicas are sometimes exchanged. The REMC method requires a tuned set of temperatures to work efficiently. The set should includes a temperature which is high enough to escape from any local minima and a temperature which is low enough to search the grand state. Additionally, a number of temperatures should be sufficient so that exchange ratios between adjacent temperatures are high enough.
Usually, a temperature set is determined by preliminary simulations and the set is fixed throughout simulations. However, many temperatures, and consequently, many replicas are required for Sudoku problems since the range of energy is wide. Therefore, we adjust the temperature set simultaneously throughout simulations to keep exchange ratios between replicas. While we cannot obtain the canonical ensemble without a fixed set of temperatures, it is not problem since we are interested only in the configuration having the lowest energy, In order to obtain a temperature set which achieves same exchange ratio between neighboring temperatures, the following procedure is proposed [10,12].
where M is a number of replicas, β n m is the inverse temperature of the m-th replica at the n-th exchange of temperature, and p n m is the acceptance ratio of exchange between the m-th and the (m + 1)-th replicas, respectively. The value of p m is estimated from MC sampling between exchange. After convergence, the acceptance ratios will share the identical value as p 1 = p 2 = · · · = p M−1 . However, with a large number of iterations, the temperature set can converge into the following trivial state In order to prevent the trivial convergence, we determine a desired value of acceptance ratiop as follows, Let N s is a number of MC steps between exchange processes, i.e., a number of samples to estimate the acceptance ratios {p n m }. If energy difference between adjacent replicas is extremely large, then the acceptance ratio between the replicas becomes extremely small. Then the acceptance ratio p n m is estimated to be zero. Once it happens, we have β m = β m+1 throughout the simulations which means that the number of replicas virtually decreases. In order to avoid the above, we adopts 1/N s instead of p n m if the exchange is not performed between m-th and (m + 1)-th replica in N s MC steps. Finally, we obtain the following procedure to adjust temperatures as, The above procedure guarantees that the all temperatures have different value and a temperature increases when some of replicas are trapped in a local minimum. We choose the initial set of temperatures as The highest temperature β 1 is fixed throughout simulations.

C. Details of Simulations
From preliminary simulations, we adopt the interaction energy J = 100 and the external field h = 1 for both depth-first and width-first order calculations. First we create a solution of Sudoku puzzle randomly, then we perform the depth-first order search with the simple MC simulation with β = 0.05. After we find a puzzle with depth larger than 8, we perform the width-first search with the REMC method. We choose a number of MC steps between temperature exchange to be N s = 100. The highest temperatures is set to be β 1 = 0.01 which is high enough to escape any local minim. The desired acceptance ratio for exchangep and a number of replicas M are chosen to bep = 0.8 and M = 10, respectively. The most time consuming part of this simulation is to calculate energy. Therefore, we adopt the normal width for internal energy in Eq. (3) instead of the average width to save computational time. We list up the candidates of hard puzzles from simulations, and determine the hardest one by calculating average width for each candidate. Since computational costs strongly depend on temperature, it is inefficient to perform parallel computation for REMC due to load imbalance. Therefore, we assign all replicas to one process, i.e., a simulation of each replica is performed serially.

IV. RESULTS
In order to investigates the efficiency of the simultaneous temperature adjustments, we perform simulations with and without the adjustments. Each simulation contains 10 replicas, and 2048 independent samples are investigated for simulations with and without adjustments. We adopt the trivial parallelization, i.e., simulations are performed independently with different seeds of random numbers. Computations are carried on SGI Altix ICE 8400EX at the Institute for Solid State Physics, the University of Tokyo. Computational time is 24 hours for each run and the total amount of the computational time is about 11 CPU-core-years. . While ten temperatures are shown, three energies out of ten are shown for the visibility. The highest temperature, i.e., the lowest inverse temperature, is fixed to be β1 = 0.01 throughout simulations.
One can see that the temperatures oscillate. Thanks to the oscillation, the replicas can escape from local minima and visit both high-and low-energy states.
The cumulative distribution function (CDF) of the minimum energies found by the simulations are shown in Fig. 3. The CDF P (E) is a probability that the lowest energy found by each process is smaller than E. It shows that the simulations with temperature adjustments is more efficient than that without adjustments.
Typical time evolutions of temperatures and energies of a run with the temperature adjustments are shown in Fig. 4. The reason that temperatures oscillate is as follows. When one replica is trapped into a local minimum, then acceptance ratio of exchange temperature between the replica and its neighbor decreases. It increases temperature of the replica. After the replica escapes from the local minimum, then the acceptance ratio increases and the temperature decreases again. From the variation rage of the lowest temperature, we estimate that over 35 replicas are necessary to keep the desired acceptance ratio without the temperature adjustments, while we used only 10 replicas. It means that the simultaneous temperature adjustments work efficiently. Figure 4 (b) shows that the energies of replicas fluctuate from highenergy states to low-energy states which means that the exchange of temperatures works efficiently.
We obtain 2048 candidates of difficult puzzles from the simulations with the simultaneous temperature adjustments. We calculate an average width for each candidates, and determine the hardest one which is shown in Fig. 5. The hardest puzzle's depth is 10 and average width is 100571 ± 1198, respectively. This is about 44 times more difficult than Inkara2012. If one adopts only pencil marks and recursive backtracking techniques to solve this puzzle, then about 50000 times recursive backtracking is necessary to solve it, and therefore, it is almost impossible to solve by hand.

V. SUMMARY AND DISCUSSION
We define the Hamiltonian describing the difficulty of Sudoku puzzles. Then creating difficult puzzles reduces to minimizing the energy defined by the Hamiltonian. We perform the REMC method to minimize the Hamiltonian, and succeed to create a Sudoku puzzle which is much harder than Inkara2012 in our definition of difficulty. To our best knowledge, this is the world hard-est puzzle ever created. While the REMC is a method to obtain canonical ensemble of different temperatures simultaneously, we propose the REMC with simultaneous adjustments of temperature which does not guarantee canonical ensembles. The results presented in this manuscript demonstrates that the REMC method is useful not only for physical problems, but for general optimization problems.
A definition of difficulty strongly depends on solving algorithms. We adopt only two techniques, pencil marks and recursive backtracking, while there are many kinds of techniques to solve Sudoku puzzles. Therefore, the difficulty defined in the present manuscript can be different from the actual feeling. But once the definition of the difficulty is given, then difficult puzzles can be created in that definition since the method proposed in the present manuscript is general and it is independent of a definition of difficulty.
While we adopt the REMC method here, there are other optimization methods, such as genetic algorithm or simulated annealing, and so forth. It is one of the further issues to compare such algorithms with REMC. Since we are interested only in the lowest energy state, it is not necessary to achieve canonical ensemble for given temperature. Therefore, we have a choice of the transition probability. While we adopt the Boltzmann weight both for MC and REMC, it is possible to accelerate finding lower energy states by adopting general transition probability [13].
In the present manuscript, we create a difficult Sudoku puzzle. It is not always true that a difficult puzzle is interesting one. However, if we can define a quantity describing how interest a puzzle is, then we can create interesting puzzles by adopting the similar method proposed here.
The programs used in the present manuscript are published as open source software [14].