当前位置：首页 > news >正文

量子退火实战（1）：用PyQUBO求解数独问题的Ising模型构建

news 2026/5/16 5:56:50

1. 数独问题与量子退火的奇妙结合

数独作为一种经典的逻辑游戏，规则简单却充满挑战。标准的9x9数独要求在每一行、每一列以及每一个3x3的小宫格内填入数字1到9且不重复。传统解法通常采用回溯算法或约束满足方法，但今天我们要尝试一种全新的思路——用量子退火来解决数独问题。

我第一次接触这个想法时也觉得不可思议，直到真正用PyQUBO实现后才明白其中的精妙。量子退火特别适合解决这类组合优化问题，因为它能高效搜索庞大的解空间。而PyQUBO这个Python库，则让我们能够用熟悉的编程方式构建Ising模型，省去了手动推导QUBO矩阵的繁琐过程。

在实际项目中，我发现将数独规则转化为Ising模型有几个关键优势：首先，约束条件可以自然地表示为二次多项式；其次，量子退火算法能够并行探索多个潜在解；最重要的是，这种方法可以推广到其他类似的约束满足问题上。下面我就详细讲解如何一步步实现这个转换过程。

2. 数独规则的数学表达

2.1 变量定义与基本约束

要用量子退火解决数独，首先需要将问题转化为适合QUBO/Ising模型的形式。我采用的变量定义方式是：为每个单元格的每个可能数字创建一个二元变量。具体来说，对于一个9x9的数独，我们需要9×9×9=729个二元变量x_{i,j,k}，其中i,j表示单元格的行列位置(1-9)，k表示可能的数字(1-9)。

这里有个实用技巧：变量可以定义为{0,1}（QUBO）或{-1,+1}（Ising）。经过多次尝试，我发现对于数独问题，使用{0,1}变量更直观，因为可以直接用1表示"该数字被选中"，0表示"未被选中"。PyQUBO支持这两种变量类型，通过vartype参数指定。

基础约束有三类：

每个单元格必须且只能选择一个数字
同一行内不能有重复数字
同一列和同一宫格内也不能有重复数字

2.2 约束条件的数学转化

将这些规则转化为数学表达式需要一些技巧。以"每个单元格必须且只能选择一个数字"为例，我们可以表示为：

sum_{k=1}^9 x_{i,j,k} = 1 对于所有i,j

在QUBO模型中，等式约束通常转化为平方差最小化的形式：

(sum_{k=1}^9 x_{i,j,k} - 1)^2

这个表达式在满足约束时值为0，否则为正数。类似地，行约束可以表示为：

sum_{i=1}^9 x_{i,j,k} = 1 对于所有j,k

转化为(sum_{i=1}^9 x_{i,j,k} - 1)^2

列约束和宫格约束的转化方式也类似。实际操作中，我发现将这些约束的权重设置得当非常重要——太弱可能导致约束被违反，太强则可能影响求解效率。

3. 使用PyQUBO构建模型

3.1 初始化问题变量

有了数学表达，现在可以用PyQUBO来实现。首先安装必要的库：

pip install pyqubo dimod

然后初始化变量数组。我习惯用字典来管理这些变量，方便后续引用：

from pyqubo import Array, Constraint # 创建9x9x9的变量数组 x = Array.create('x', shape=(9,9,9), vartype='BINARY')

3.2 构建目标函数

接下来将各种约束组合成目标函数。这里有个实用建议：为不同类型的约束设置不同的权重系数，这样可以在求解时更好地平衡约束满足和优化目标。

# 单元格约束：每个格子必须且只能选一个数字 cell_constraints = sum( Constraint((sum(x[i,j,k] for k in range(9)) - 1)**2, label=f'cell_{i}_{j}') for i in range(9) for j in range(9) ) # 行约束：每行数字不重复 row_constraints = sum( Constraint((sum(x[i,j,k] for i in range(9)) - 1)**2, label=f'row_{j}_{k}') for j in range(9) for k in range(9) ) # 列约束：每列数字不重复 col_constraints = sum( Constraint((sum(x[i,j,k] for j in range(9)) - 1)**2, label=f'col_{i}_{k}') for i in range(9) for k in range(9) ) # 宫格约束：每个3x3小宫格内数字不重复 box_constraints = sum( Constraint((sum(x[i+3*bi, j+3*bj, k] for i in range(3) for j in range(3)) - 1)**2, label=f'box_{bi}_{bj}_{k}') for bi in range(3) for bj in range(3) for k in range(9) ) # 组合所有约束，设置适当权重 H = 1.0*cell_constraints + 1.0*row_constraints + 1.0*col_constraints + 1.0*box_constraints

3.3 处理已知数字

实际数独问题通常有部分已知数字。这些可以作为固定约束直接加入模型：

# 假设已知(0,0)位置数字为5，(4,4)位置数字为3 known_values = {(0,0):5, (4,4):3} # 添加已知值约束 for (i,j), k in known_values.items(): H += Constraint((1 - x[i,j,k-1])**2, label=f'known_{i}_{j}') * 2.0

这里我给了已知值约束更高的权重(2.0)，确保它们被优先满足。

4. 模型求解与结果解析

4.1 编译模型并采样

现在可以将模型编译为BQM(二元二次模型)并求解了：

model = H.compile() bqm = model.to_bqm() # 使用模拟退火求解器 from neal import SimulatedAnnealingSampler sampler = SimulatedAnnealingSampler() sampleset = sampler.sample(bqm, num_reads=1000)

在实际测试中，我发现num_reads参数对结果质量影响很大。对于9x9数独，通常需要1000次以上的读取才能稳定找到有效解。

4.2 解码和验证结果

得到采样结果后，需要解码并验证：

# 解码样本 decoded_samples = model.decode_sampleset(sampleset) # 找到能量最低(最可能满足约束)的样本 best_sample = min(decoded_samples, key=lambda x: x.energy) # 检查约束是否满足 print("违反的约束:", best_sample.constraints(only_broken=True))

如果发现有约束被违反，可能需要调整权重或增加采样次数。理想情况下，我们应该得到一个能量为0且没有约束违反的解。

4.3 可视化数独解

最后，将解转换为更友好的数独格式：

import numpy as np # 初始化9x9数独网格 sudoku_grid = np.zeros((9,9), dtype=int) # 填充解 for i in range(9): for j in range(9): for k in range(9): if best_sample.array('x', (i,j,k)) == 1: sudoku_grid[i,j] = k+1 print("解得的数独:") print(sudoku_grid)