数独を疑似量子アニーリングで解く

はじめに

量子アニーリング（Quantum Annealing）やシミュレーテッドアニーリング（Simulated Annealing）は、組合せ最適化問題を解く手法としてよく登場する。
Pythonライブラリ OpenJij を使って数独をQUBOに変換し、疑似量子アニーリングで解けるか試してみた。

結論: 簡単な盤面は解けたが、「世界一難しい数独」と呼ばれる問題は解けなかった。

実行環境

• OS: Ubuntu
• 言語：Python 3.11

使用ライブラリ

• dimod
  QUBO（Quadratic Unconstrained Binary Optimization）モデルを定義するライブラリ。D-Wave関連でよく使われる。

• OpenJij
  疑似量子アニーリングやシミュレーテッドアニーリングをPythonから実行できるライブラリ。

  • SASampler: 古典的シミュレーテッドアニーリング
  • SQASampler: 疑似量子アニーリング

QUBOモデル化

今回ベースとなるコードは生成AIで出力し、結果を確認しやすいように改良を加えた。

数独のルールをQUBOに落とし込む方法:

• セル制約: 各マスは1〜9のどれか1つだけ
• 行制約 / 列制約: 各行・各列で同じ数字は1回だけ
• ブロック制約: 3x3ブロック内で同じ数字は1回だけ
• ヒント: 与えられた数字を強めのバイアスで固定する

$\sum_{} x=1$ を理想として、これらを ($1 - \sum_{}x^2$) の形でペナルティ項にして、BQM（Binary Quadratic Model）に変換する。

実装コード例

solved, energy, sampleset = solve_sudoku_with_openjij(
    puzzle,
    num_reads=100,
    sampler_type="sa",  # "sqa" に変えると疑似量子アニーリング
    A=4.0,             # 制約の重み
    C=10.0              # ヒントの重み
)

• num_reads: サンプリング回数
• sampler_type: "sa" or "sqa"
• A, C: 制約とヒントの強さ

基本、num_readsとA, Cの数値を増加させてsampler_typeを"sqa"にすると精度が上がるようだ。

結果の確認

結果が出力されたら数独のルールに則っているかをチェックするコードを実行し、整合性を確認する。

def check_sudoku_solution(grid):
    """
    grid: 9x9 数独盤
    ルール違反があればその位置と内容を返す
    正しければ True を返す
    """
    errors = []

    # 行チェック
    for r in range(9):
        seen = set()
        for c in range(9):
            v = grid[r][c]
            if v < 1 or v > 9:
                errors.append(f"行 {r+1} 列 {c+1}: 値 {v} が1-9の範囲外")
            elif v in seen:
                errors.append(f"行 {r+1} に重複: {v}")
            else:
                seen.add(v)

    # 列チェック
    for c in range(9):
        seen = set()
        for r in range(9):
            v = grid[r][c]
            if v in seen:
                errors.append(f"列 {c+1} に重複: {v}")
            else:
                seen.add(v)

    # 3x3 ブロックチェック
    for br in range(3):
        for bc in range(3):
            seen = set()
            for dr in range(3):
                for dc in range(3):
                    r, c = br*3 + dr, bc*3 + dc
                    v = grid[r][c]
                    if v in seen:
                        errors.append(f"ブロック ({br+1},{bc+1}) に重複: {v}")
                    else:
                        seen.add(v)

    if errors:
        print("違反が見つかりました:")
        for e in errors:
            print(" -", e)
        return False
    else:
        print("正しい解答です！")
        return True

コード全体(Python環境があれば動かしてみてください)

import itertools
from collections import defaultdict
import dimod
import openjij as oj
from tqdm import tqdm


def build_sudoku_bqm(grid, A=4.0, C=10.0):
    """
    grid: 9x9 list of lists, 0 = blank, 1..9 = given numbers
    A: penalty for Sudoku constraints
    C: bias strength for clues (given cells)
    """

    # 変数キーを文字列で統一
    def key(r, c, n):
        # r,c,n are 0-indexed; n in 0..8
        return f"x_{r}_{c}_{n}"

    # QUBO（上三角も含めて辞書で管理）
    Q = defaultdict(float)
    offset = 0.0

    # 汎用：one-hot  (1 - sum_i x_i)^2
    def add_one_hot(keys, weight):
        nonlocal offset
        offset += weight  # 1 * weight
        # 線形項
        for i in keys:
            Q[(i, i)] += -1.0 * weight
        # 二次項（i<j）
        for i, j in itertools.combinations(keys, 2):
            Q[(i, j)] += 2.0 * weight

    # --- 変数の用意（9x9x9=729） ---
    all_keys = []
    for r in range(9):
        for c in range(9):
            for n in range(9):
                all_keys.append(key(r, c, n))

    # --- セルごとに one-hot：各マスは 1..9 のいずれか1つ ---
    for r in range(9):
        for c in range(9):
            add_one_hot([key(r, c, n) for n in range(9)], A)

    # --- 行制約：各行の各数字は1回だけ ---
    for r in range(9):
        for n in range(9):
            add_one_hot([key(r, c, n) for c in range(9)], A)

    # --- 列制約：各列の各数字は1回だけ ---
    for c in range(9):
        for n in range(9):
            add_one_hot([key(r, c, n) for r in range(9)], A)

    # --- ブロック制約：3x3 ブロック内で各数字は1回だけ ---
    for br in range(3):
        for bc in range(3):
            cells = [(br * 3 + dr, bc * 3 + dc) for dr in range(3) for dc in range(3)]
            for n in range(9):
                add_one_hot([key(r, c, n) for (r, c) in cells], A)

    # --- 与えられたヒントを強く“推す”バイアス（x_{rcn}=1 を好む）---
    for r in range(9):
        for c in range(9):
            v = grid[r][c]
            if v != 0:
                # v は 1..9、n は 0..8 なので n = v-1
                Q[(key(r, c, v - 1), key(r, c, v - 1))] += -C

    # dimod の BQM を作成
    bqm = dimod.BinaryQuadraticModel.from_qubo(dict(Q), offset=offset)
    return bqm


def solve_sudoku_with_openjij(grid, num_reads=200, sampler_type="sqa", A=4.0, C=10.0):
    batch_size = int(num_reads / 10)
    """
    grid: 数独盤 (9x9, 0=空白)
    num_reads: サンプル数（例: 200）
    batch_size: 一度にサンプルする回数（例: 50）
    sampler_type: "sa" or "sqa"
    """
    bqm = build_sudoku_bqm(grid, A, C)

    if sampler_type == "sqa":
        sampler = oj.SQASampler()
    else:
        sampler = oj.SASampler()

    # サンプルを分割して進捗バーを出す
    all_samples = None
    for i in tqdm(range(0, num_reads, batch_size), desc="Annealing Progress"):
        current_batch = min(batch_size, num_reads - i)
        sampleset = sampler.sample(bqm, num_reads=current_batch)
        if all_samples is None:
            all_samples = sampleset
        else:
            all_samples = dimod.concatenate([all_samples, sampleset])

    # 最良解を取り出す
    best = all_samples.first.sample
    solved = [[0] * 9 for _ in range(9)]
    for r in range(9):
        for c in range(9):
            # 元のヒントがある場合はそのまま
            if grid[r][c] != 0:
                solved[r][c] = grid[r][c]
                continue

            chosen_n = None
            best_val = -1
            for n in range(9):
                v = best.get(f"x_{r}_{c}_{n}", 0)
                if v > best_val:
                    best_val = v
                    chosen_n = n
            solved[r][c] = chosen_n + 1  # 1..9 に戻す

    return solved, all_samples.first.energy, all_samples


def print_grid(g):
    for r in range(9):
        row = ""
        for c in range(9):
            row += f"{g[r][c]} "
            if c % 3 == 2 and c != 8:
                row += "| "
        print(row.strip())
        if r % 3 == 2 and r != 8:
            print("-" * 21)


def check_sudoku_solution(grid):
    """
    grid: 9x9 数独盤
    ルール違反があればその位置と内容を返す
    正しければ True を返す
    """
    errors = []

    # 行チェック
    for r in range(9):
        seen = set()
        for c in range(9):
            v = grid[r][c]
            if v < 1 or v > 9:
                errors.append(f"行 {r + 1} 列 {c + 1}: 値 {v} が1-9の範囲外")
            elif v in seen:
                errors.append(f"行 {r + 1} に重複: {v}")
            else:
                seen.add(v)

    # 列チェック
    for c in range(9):
        seen = set()
        for r in range(9):
            v = grid[r][c]
            if v in seen:
                errors.append(f"列 {c + 1} に重複: {v}")
            else:
                seen.add(v)

    # 3x3 ブロックチェック
    for br in range(3):
        for bc in range(3):
            seen = set()
            for dr in range(3):
                for dc in range(3):
                    r, c = br * 3 + dr, bc * 3 + dc
                    v = grid[r][c]
                    if v in seen:
                        errors.append(f"ブロック ({br + 1},{bc + 1}) に重複: {v}")
                    else:
                        seen.add(v)

    if errors:
        print("違反が見つかりました:")
        for e in errors:
            print(" -", e)
        return False
    else:
        print("正しい解答です！")
        return True


if __name__ == "__main__":
    # 例：完成盤から穴を空けた問題（0 が空白）好きな問題を入れる
    puzzle = [
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0],
    ]

    print("Puzzle:")
    print_grid(puzzle)

    solved, energy, sampleset = solve_sudoku_with_openjij(puzzle, num_reads=100, sampler_type="sa", A=4.0, C=8.0)
    print("\nSolved (energy =", energy, "):")
    print_grid(solved)
    is_correct = check_sudoku_solution(solved)

実行結果

簡単な盤面

簡単と言っても、生成AIに「最難関の数独を作って」と命令して出力されたパズル。

    puzzle = [
        [0, 0, 0, 0, 0, 0, 0, 3, 0],
        [0, 0, 0, 0, 0, 7, 0, 0, 0],
        [0, 0, 0, 6, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, 0, 0],
        [0, 7, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 4, 0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 0, 0, 2],
        [0, 0, 0, 0, 8, 0, 0, 0, 0],
        [1, 0, 0, 0, 0, 0, 0, 0, 0],
    ]

このような問題はsaでも十分解けた。(穴が空きすぎて逆に解の自由度が高そう)

Solved (energy = -90.0 ):
9 2 7 | 5 4 1 | 6 3 8
6 3 1 | 8 2 7 | 4 9 5
8 4 5 | 6 3 9 | 7 2 1
---------------------
5 6 8 | 9 7 2 | 1 4 3
4 7 9 | 3 1 5 | 2 8 6
2 1 3 | 4 6 8 | 9 5 7
---------------------
3 8 4 | 1 9 6 | 5 7 2
7 5 6 | 2 8 4 | 3 1 9
1 9 2 | 7 5 3 | 8 6 4
正しい解答です！

制約違反もなく正しい解答になった。
以下のパラメータで10秒以内に解が出る。※間違えることもある

• num_reads: 100
• sampler_type: "sa"
• A=4,C=10

世界一難しいと呼ばれる数独

「世界一難しい数独」と呼ばれる盤面を試したが解けなかった。
①https://rocketnews24.com/2012/07/03/22654/

    puzzle = [
    [8, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 3, 6, 0, 0, 0, 0, 0],
    [0, 7, 0, 0, 9, 0, 2, 0, 0],
    [0, 5, 0, 0, 0, 7, 0, 0, 0],
    [0, 0, 0, 0, 4, 5, 7, 0, 0],
    [0, 0, 0, 1, 0, 0, 0, 3, 0],
    [0, 0, 1, 0, 0, 0, 0, 6, 8],
    [0, 0, 8, 5, 0, 0, 0, 1, 0],
    [0, 9, 0, 0, 0, 0, 4, 0, 0],
    ]

②https://gigazine.net/news/20100822_hardest_sudoku/

    puzzle = [
    [0, 0, 5, 3, 0, 0, 0, 0, 0],
    [8, 0, 0, 0, 0, 0, 0, 2, 0],
    [0, 7, 0, 0, 1, 0, 5, 0, 0],
    [4, 0, 0, 0, 0, 5, 3, 0, 0],
    [0, 1, 0, 0, 7, 0, 0, 0, 6],
    [0, 0, 3, 2, 0, 0, 0, 8, 0],
    [0, 6, 0, 5, 0, 0, 0, 0, 9],
    [0, 0, 4, 0, 0, 0, 0, 3, 0],
    [0, 0, 0, 0, 0, 9, 7, 0, 0],
    ]

▼重複が残り、制約違反が出た。

Solved (energy = -176.0 ):
2 4 5 | 3 8 6 | 1 9 7
8 3 1 | 9 5 7 | 6 2 4
6 7 9 | 4 1 2 | 5 3 8
---------------------
4 8 6 | 1 9 5 | 3 7 2
9 1 2 | 8 7 3 | 4 5 6
7 5 3 | 2 6 4 | 9 8 1
---------------------
1 6 7 | 5 3 8 | 2 4 9
5 9 4 | 7 2 1 | 8 3 3
3 2 8 | 6 4 9 | 7 1 5
違反が見つかりました:
 - 行 8 に重複: 3
 - 列 8 に重複: 3
 - ブロック (3,3) に重複: 3

色々と試行錯誤し、パラメータを以下のように調整して30分くらい処理時間を要したが、それでも解けなかった。

• num_reads: 5000
• sampler_type: "sqa"
• A=16,C=30

所感

• 簡単な問題は解ける
  → 制約が弱く収束しやすい(正解となるパターンが多い)
• 難しい問題(制約が多い問題)は解きにくい
  → 難しい問題にはちゃんと難しい理由があるのだと感じた

余談だが、「sudoku solver」などで検索すると出るサイトで上記世界一難しい数独を入れると一瞬で解けるので、解法が確立している問題については古典的アルゴリズムの方が早く解けてしまう問題もあるのだと実感した。

今後やってみたいこと

• パラメータを変えたときの傾向を知りたい
• 量子コンピュータ実機で試す
• 数独以外の組合せ最適化問題も試す