Skip to main content
Have a personal or library account? Click to login
A Quantum-Inspired Algorithm for Solving Sudoku Puzzles and the MaxCut Problem † Cover

A Quantum-Inspired Algorithm for Solving Sudoku Puzzles and the MaxCut Problem †

By:  and    
Open Access
|Dec 2025

Figures & Tables

Figure 1.

The spin representations of three Sudoku puzzles published in The New York Times on January 2 (diamond), January 12 (circle), and January 14 (plus), 2025, are shown. The mean and standard deviation of the magnetic field hmz are as follows: 8.7 ± 3.0 for the January 2 puzzle (diamond), 8.0 ± 3.7 for January 12 (circle), and 7.3 ± 2.7 for January 14 (plus).

Figure 2.

The dense, long-range couplings among the 210 spins for the January 12, 2025, Sudoku puzzle. Each line represents an Ising coupling Jm,n = 2, which favors spin-m and spin-n being antiparallel. Each spin can point either up or down. Finding the lowest-energy spin configuration corresponds to solving the Sudoku puzzle.

Figure 3.

Statistical analysis of the couplings between the spins. For nearest neighbors, about 75% are coupled. The coupling fraction drops to the 10% level when the distance grows to 6 or 7, then fluctuates between 1% and 10%. The legends are the same as Figure 1. The vertical axis is in logarithmic scale.

Table 1.

A few examples of solved Sudoku puzzles in chronological order according to their publication date in The New York Times. These puzzles vary in the number of clues (Nc), total number of spins (Ns), and other parameters; see the main text for details.

Puzzle dateLevelNcNNsc1+c2+c3
January 8, 2024h2457211368.5
December 10, 2024m2655200359.0
January 2, 2025m2259250520.5
January 8, 2025m2358232461.5
January 12, 2025m2754210426.5
January 14, 2025m2358232457.0
January 15, 2025m2556210380.0
Figure 4.

How the January 8, 2024, puzzle gets solved: As the driving field hx is reduced, the energy gradually drops to zero. The total spin along the x-direction vanishes, while the total spin along the z-direction approaches its ground state value. The bond dimension is D = 20, and the driving field is set to hx = 0.75.

Figure 5.

The solution process for another puzzle with fewer clues and more spins: D = 20, hx = 1.

Figure 6.

The nonlinear drop in energy during the driving process. See main text for details.

Figure 7.

The evolution of the 251 spins during the solution process of the January 2, 2025, puzzle. The value of Szm is color-coded: red represents spin-up, and blue represents spin-down.

281354679
973816524
564792813
139468257
825973146
647521938
492635781
356187492
718249365
Table 2.

Twenty weighted MaxCut instances and the parameters used to solve them. Each instance is identified by its file name in the Biq Mac Library. The edge weights are either 1 or – 1. For the first ten instances, the number of vertices NV = 80 and the number of edges NE = 316. For the next ten instances, NV = 100 and NE = 495. The algorithm successfully finds E0, the ground state energy of the spin model (the MaxCut value is simply – E0), except for the case pm1s_100.5.

InstanceE0MhxηD
pm1s_80.0–7951030
pm1s_80.1–6951030
pm1s_80.2–67510.330
pm1s_80.3–6651030
pm1s_80.4–69510.330
pm1s_80.5–662010.330
pm1s_80.6–7151030
pm1s_80.7–6951030
pm1s_80.8–6851030
pm1s_80.9–6751030
pm1s_100.0–127510.330
pm1s_100.1–126510.330
pm1s_100.2–1251010.330
pm1s_100.3–1111010.330
pm1s_100.4–128510.330
pm1s_100.5125101060
pm1s_100.6–1221010.330
pm1s_100.7–1122010.360
pm1s_100.8–1201010.330
pm1s_100.9–127510.330
Figure 8.

(a) Graph representation of the problem. The edges carry weights of 1 (in green) or – 1 (in blue). The solution is indicated by the node color: nodes in black belong to set A, while the rest of the nodes (in red) belong to set Ā. (b) The energy (dot), expectation values of Sz (cross, measured from its ground state value and magnified by 10 times), and Sx (circle) during the driving process. Solving the weighted MaxCut instance pm1s_100.9 from the Biq Mac Library.

Table 3.

Twenty unweighted MaxCut instances from the Biq Mac Library with dense edges. The first ten cases have NV = 60 and NE = 885, while the remaining cases have NV = 100 and NE = 2,475. All edge weights are set to 1. The algorithm successfully solves all instances except g05_100.3 and g05_100.5, where the E0 values should be – 1,424 and – 1,436, respectively [39,40]. Here, Nsw refers to the number of DMRG sweeps per driving step, with hx = 1 and η = 0.3.

InstanceE0MNswD
g05_60.0–53610530
g05_60.1–53210540
g05_60.2–52910530
g05_60.3–53810530
g05_60.4–52710530
g05_60.5–53310530
g05_60.6–53110530
g05_60.7–53510560
g05_60.8–53010530
g05_60.9–53320540
g05_100.0–143020540
g05_100.1–142520540
g05_100.2–143220540
g05_100.3–142310530
g05_100.4–144020540
g05_100.5–143510530
g05_100.6–1434101040
g05_100.7–1431101040
g05_100.8–1432101060
g05 100.9–1430101040
Figure 9.

(a) Graph representation of the problem. All edges (green lines) carry equal weight of 1. The MaxCut found by our algorithm is color-coded: nodes in black belong to set Ā, while nodes in red belong to set Ā. (b) The expectation value of Szm during the driving process as the spins gradually settle into their ground state orientation (red represents spin up, blue represents spin down). Only driving steps 2 to 10 are shown to enhance the color contrast. Solving the unweighted MaxCut instance g05_100.4 from the Biq Mac library.

Table 4.

The algorithm parameters used to successfully solve five MaxCut instances from the Biq Mac Library with 251 vertices and over 3200 edges, each with integer weights; D = 30, M = 10.

InstanceE0NVNEhxηNsw
bqp250-2-4481025132850.05010
bqp250-4-4127425133970.05010
bqp250-6-4101425134330.3005
bqp250-8-3572625132650.05010
bqp250-10-4044225132941.000.35
Figure 10.

(a) The distribution of the edge weights wi,j sorted by |ij|, the distance between the nodes. Most of the weights lie within the interval from –100 to 100. (b) The variation of energy during the driving process for three different values of hx. All three choices converge to the ground state, yielding the correct MaxCut value. The edge matrix has been normalized such that the maximum absolute value of its elements is 1. Solving MaxCut instance bqp250-8 with 251 nodes and 3265 edges.

repeat
 hop, sweep
until M times
DOI: https://doi.org/10.2478/qic-2025-0023 | Journal eISSN: 3106-0544 (formerly 1533-7146) | Journal ISSN: 1533-7146
Language: English
Page range: 397 - 420
Submitted on: Jun 5, 2025
Accepted on: Aug 6, 2025
Published on: Dec 31, 2025
Published by: Cerebration Science Publishing Co., Limited
In partnership with: Paradigm Publishing Services
Publication frequency: 1 issue per year

© 2025 Max B. Zhao, Fei Li, published by Cerebration Science Publishing Co., Limited
This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.