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: Max B. Zhao and  Fei Li  
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 hmzh_m^z 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).
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 hmzh_m^z 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.
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.
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.

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.
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.
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.
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 SzmS_m^z is color-coded: red represents spin-up, and blue represents spin-down.
The evolution of the 251 spins during the solution process of the January 2, 2025, puzzle. The value of SzmS_m^z is color-coded: red represents spin-up, and blue represents spin-down.

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.
(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.

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 SzmS_z^m 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.
(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 SzmS_z^m 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.

Figure 10.

(a) The distribution of the edge weights wi,j sorted by |i – j|, 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.
(a) The distribution of the edge weights wi,j sorted by |i – j|, 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.

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

j_qic-2025-0023_utab_002

repeat
 hop, sweep
until M times

j_qic-2025-0023_utab_001

281354679
973816524
564792813
139468257
825973146
647521938
492635781
356187492
718249365

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

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

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
DOI: https://doi.org/10.2478/qic-2025-0023 | Journal eISSN: 3106-0544 | 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
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.