第13章 QUBOで数独を解く

数独の規則を、守るか破るかで点数が変わる式にしてみます。正しい完成盤面なら0点、同じ行に 同じ数字を置くなどの違反があれば加点される式です。あとは、この点数を最小にする0と1の組合せを 探します。

このようなモデルを QUBO（Quadratic Unconstrained Binary Optimization） と呼びます。日本語にすれば「制約なし二値二次最適化」です。QUBOは問題の書き方であり、 解法や計算機を指定しません。掲載コードでは、通常のCPU上で焼きなまし法を実行します。

一つだけ選ぶ条件を点数にする

まず、二つの候補 \(x_1,x_2\) から一つだけ選ぶ場面を考えます。選ぶ候補を1、選ばない候補を 0とすれば、次の式を罰点にできます。

\[P(x_1,x_2)=(x_1+x_2-1)^2\]

一つだけ選んだときは括弧の中が0になります。両方とも選ばない場合と、両方とも選ぶ場合は1です。 四通りをプログラムで確かめると、次の表になります。

"""「二変数から一つだけ選ぶ」罰則の全入力を表示します。"""


def penalty(first: int, second: int) -> int:
    """Return ``(first + second - 1)^2``."""
    return (first + second - 1) ** 2


def main() -> None:
    print("x1 x2 penalty")
    for first in (0, 1):
        for second in (0, 1):
            # ちょうど一方が1のときだけ、違反を表す罰点が0になる。
            print(first, second, penalty(first, second))


if __name__ == "__main__":
    main()

x1 x2 penalty
0 0 1
0 1 0
1 0 0
1 1 1

この式には二次式より高い項がありません。さらに、0か1を取る変数では \(x_i^2=x_i\) なので、 展開すると次のQUBOになります。

\[P(x_1,x_2)=1-x_1-x_2+2x_1x_2\]

一般に、候補 \(x_1,\ldots,x_n\) から一つだけ選ぶ条件も、同じ形で書けます。

\[\left(\sum_{i=1}^n x_i-1\right)^2 =1-\sum_{i=1}^n x_i+2\sum_{i<j}x_ix_j\]

右辺には、各変数だけに掛かる一次の項と、二変数の組に掛かる二次の項しかありません。条件を ソルバーへ別に渡す代わりに、破ったときに大きくなる式へ埋め込んでいます。このため unconstrained と呼ばれます。ただし、変数が0か1である制限は残っています。

QUBOは何を最小化するか

QUBOでは、二値変数のベクトル \(\boldsymbol{x}\) に対し、一次と二次の項を足したエネルギー \(E(\boldsymbol{x})\) を最小化します。ここでいうエネルギーは、候補の良し悪しを比べる点数です。 低いほど罰点が少なくなります。この章では係数の意味が見えるよう、次の書き方をします。

\[E(\boldsymbol{x}) =C+\sum_i a_i x_i+\sum_{i<j}b_{ij}x_ix_j, \qquad x_i\in\{0,1\}\]

\(C\) は定数、\(a_i\) は一変数の係数、\(b_{ij}\) は二変数の組の係数です。文献や ライブラリでは行列 \(Q\) を使って \(\boldsymbol{x}^{\mathsf T}Q\boldsymbol{x}\) と書く こともあります。対称な成分をどう数えるかで係数の置き方が変わるため、コードでは dimod.BinaryQuadraticModel の一次係数、二次係数、定数項を分けて渡します。

QUBOは、割当て、巡回路、スケジューリングなどの組合せ最適化を二値変数と罰則へ変換するときにも 使われます。制約を二次式へ移す手順と多くの定式化例は、Gloverらのチュートリアルにまとめられて います [1]。

QUBOは、古典的な厳密ソルバー、ヒューリスティック、量子アニーリングのいずれでも扱えます。 モデルを作った後、最低エネルギーを見つける処理は選んだソルバーが担当します。

数独を729個の二値変数で表す

9×9数独の行 \(r\)、列 \(c\) に数字 \(d\) を入れる候補を \(x_{r,c,d}\) とします。各マスに九候補があるので、変数は \(9\times9\times9=729\) 個です。次の四種類の条件を罰則にします。

• 一つのマスには、数字が一つだけ入ります。

• 一つの行には、各数字が一度だけ現れます。

• 一つの列には、各数字が一度だけ現れます。

• 一つの3×3ブロックには、各数字が一度だけ現れます。

各条件は「九候補から一つだけ選ぶ」ので、先ほどの二次罰則をそのまま使えます。コードでは add_exactly_one が、一次係数へ \(-1\)、候補ペアの二次係数へ \(2\)、定数項へ \(1\) を加えます。weight は罰則全体に掛ける正の重みです。

たとえば左上のマスには、 \((x_{1,1,1}+x_{1,1,2}+\cdots+x_{1,1,9}-1)^2\) を加えます。 add_exactly_one は、この式を展開した結果を linear、quadratic、定数項の三部分へ分けて 登録します。式を別の条件へ変えたのではなく、ソルバーが受け取る係数表へ移しています。

def add_exactly_one(
    linear: dict[Variable, float],
    quadratic: dict[tuple[Variable, Variable], float],
    choices: Iterable[Variable],
    weight: float,
) -> float:
    """Add ``weight * (sum(choices) - 1)^2`` and return its offset."""
    variables = tuple(choices)

    # x^2 = x なので、各候補の一次係数は -weight になる。
    for item in variables:
        linear[item] += -weight

    # 二候補を同時に選んだ分は、2 * weight * x_i * x_j で数える。
    for first, second in combinations(variables, 2):
        pair = tuple(sorted((first, second)))
        quadratic[pair] += 2 * weight

    # 展開前の式に含まれる定数項 weight を返す。
    return weight

初期配置も罰則にします。指定された候補を \(x=1\) にする罰則は \((x-1)^2\) です。二値変数なら \(1-x\) になるので、定数項に1、対応する一次係数に \(-1\) を加えます。

四種類の規則と初期配置の罰則をすべて足したものが、数独のエネルギーです。

def build_qubo(givens: Board, weight: float = 1.0) -> dimod.BinaryQuadraticModel:
    if len(givens) != SIZE * SIZE:
        raise ValueError("この例は9×9数独だけを扱います")

    linear: defaultdict[Variable, float] = defaultdict(float)
    quadratic: defaultdict[tuple[Variable, Variable], float] = defaultdict(float)
    offset = 0.0

    # 各マスでは、九つの数字候補から一つだけを選ぶ。
    for row in range(SIZE):
        for column in range(SIZE):
            offset += add_exactly_one(
                linear,
                quadratic,
                (variable(row, column, digit) for digit in range(1, SIZE + 1)),
                weight,
            )

    # 各行では、それぞれの数字を一度だけ選ぶ。
    for row in range(SIZE):
        for digit in range(1, SIZE + 1):
            offset += add_exactly_one(
                linear,
                quadratic,
                (variable(row, column, digit) for column in range(SIZE)),
                weight,
            )

    # 各列でも、それぞれの数字を一度だけ選ぶ。
    for column in range(SIZE):
        for digit in range(1, SIZE + 1):
            offset += add_exactly_one(
                linear,
                quadratic,
                (variable(row, column, digit) for row in range(SIZE)),
                weight,
            )

    # 各3×3ブロックでも、それぞれの数字を一度だけ選ぶ。
    for top in range(0, SIZE, BOX):
        for left in range(0, SIZE, BOX):
            for digit in range(1, SIZE + 1):
                offset += add_exactly_one(
                    linear,
                    quadratic,
                    (
                        variable(row, column, digit)
                        for row in range(top, top + BOX)
                        for column in range(left, left + BOX)
                    ),
                    weight,
                )

    # 初期配置の候補は1にする。(x - 1)^2 は二値変数なら 1 - x になる。
    for index, digit in enumerate(givens):
        if digit:
            row, column = divmod(index, SIZE)
            linear[variable(row, column, digit)] += -weight
            offset += weight

    return dimod.BinaryQuadraticModel(linear, quadratic, offset, dimod.BINARY)

今回は数独の規則以外に最小化したい目的がないため、すべての重みを1にしています。どの罰則も 平方で0以上なので、合計が0なら全条件を満たします。逆に、全条件を満たす盤面のエネルギーは0です。 したがって、ゼロエネルギー状態と完成盤面が対応します。

配送距離など別の目的も同時に最小化するQUBOでは、制約を破って目的値だけを小さくする状態に負けない よう、罰則の重みを決める必要があります。数独ではその競合がないため、重み選びの問題を避けられます。

初期配置で決まる候補を固定する

729変数のモデルには、同じ変数ペアをまとめると10,206個の二次項があります。 このまますべてのビットを動かすと、初期配置に反する状態まで探索します。そこで、焼きなましに 渡す前に、初期配置だけで値が決まる変数を固定します。

たとえば左上に5が指定されていれば、\(x_{1,1,5}=1\)、同じマスのほかの候補は0です。 さらに、同じ行・列・ブロックの空きマスには5を置けないので、それらの候補も0にします。 ここでは初期配置から直接分かる固定値だけを使い、空きマスの候補を順に確定する探索は行いません。

def fixed_candidates(givens: Board) -> dict[Variable, int]:
    fixed: dict[Variable, int] = {}
    for row in range(SIZE):
        for column in range(SIZE):
            given = givens[row * SIZE + column]
            top, left = row // BOX * BOX, column // BOX * BOX
            used = set(givens[row * SIZE : (row + 1) * SIZE])
            used.update(givens[column::SIZE])
            used.update(
                givens[r * SIZE + c]
                for r in range(top, top + BOX)
                for c in range(left, left + BOX)
            )
            for digit in range(1, SIZE + 1):
                if given:
                    # 初期配置のマスは、指定数字だけを1に固定する。
                    fixed[variable(row, column, digit)] = int(digit == given)
                elif digit in used:
                    # 同じ行・列・ブロックの初期配置と衝突する候補は0にする。
                    fixed[variable(row, column, digit)] = 0
    return fixed

bqm.fix_variables は固定値を式へ代入し、係数と定数項をまとめ直します。 たとえば \(2xy\) で \(x=1\) なら \(2y\)、\(x=0\) なら0になります。 固定した範囲では代入前後のエネルギーが等しく、数独の解を除くこともありません。 通常問題では、焼きなましが動かす変数は223個です。出力では、モデル全体の変数数を variables、固定後に残った変数数を sampled_variables として表示します。

焼きなましで状態をサンプリングする

モデルの構築にはD-Wave Oceanの dimod、探索には nealのSimulatedAnnealingSampler を使います。nealは、温度に相当する値を下げながらビットを反転し、低いエネルギーへ移る古典的な 焼きなまし法です。温度が高いうちは、エネルギーが上がる移動もときどき受け入れ、同じ場所に 留まり続けるのを避けます。温度が下がるにつれて、そのような移動を受け入れにくくします。 量子ハードウェアやクラウド接続は使いません。

一回の read は、焼きなましを一度実行して一状態を得る処理です。sweep では全変数に対する 更新を一巡します。nealの公式文書でも、num_reads は実行回数、num_sweeps は各実行の 更新巡回数として定義されています [2]。

この章では乱数シードを 20260808、読み出しを500回、各読み出しを2,000スイープに固定しました。 得られた状態に固定済みの候補を戻し、盤面へ復元します。エネルギーが0で、さらに共通検証器が 受理した盤面だけを解として残します。

def sample_sudoku(givens: Board, *, seed: int, reads: int, sweeps: int) -> SamplingResult:
    bqm = build_qubo(givens)
    fixed = fixed_candidates(givens)
    # 固定値を代入し、一次・二次係数と定数項をまとめ直す。
    bqm.fix_variables(fixed)
    sampler = neal.SimulatedAnnealingSampler()

    # 同じ乱数シード、読み出し回数、スイープ数なら実行条件を再現できる。
    samples = sampler.sample(
        bqm,
        seed=seed,
        num_reads=reads,
        num_sweeps=sweeps,
    ).aggregate()

    valid_boards: set[Board] = set()
    zero_energy_reads = 0
    for datum in samples.data(fields=["sample", "energy", "num_occurrences"]):
        if abs(datum.energy) > 1e-9:
            continue
        zero_energy_reads += datum.num_occurrences

        sample = dict(datum.sample)
        sample.update(fixed)
        board = decode_sample(sample)
        if board is None:
            continue
        valid, _ = validate_solution(board, givens)
        if valid:
            valid_boards.add(board)

    # 焼きなましでゼロエネルギー解を見つけられなくても、unsatとは断定しない。
    status = "solved" if valid_boards else "unknown"
    return SamplingResult(
        status=status,
        best_energy=float(samples.first.energy),
        zero_energy_reads=zero_energy_reads,
        solutions=tuple(sorted(valid_boards)),
        sampled_variables=len(bqm.variables),
    )

焼きなまし法は、限られた反復回数で最小エネルギーへ到達するとは限りません。 エネルギー0の盤面を得られなかった場合は、解なしと断定せず unknown を返します。 シードを固定して実行条件をそろえても、最適性の証明にはなりません。

実行と検証

モデルの実行および検証手順については examples/12-qubo/README.md を参照してください。

通常問題

入力は fixtures/standard-9x9.sdk です。リポジトリ直下で次を実行します。

500600902
000105308
000000500
800001020
000003000
010920800
060500004
280000000
305000070

$ uv run --frozen python examples/12-qubo/solve.py fixtures/standard-9x9.sdk \
    --seed 20260808 --reads 500 --sweeps 2000 --limit 2

status: solved
sampler: neal.SimulatedAnnealingSampler
variables: 729
sampled_variables: 223
seed: 20260808
reads: 500
sweeps: 2000
best_energy: 0.0
zero_energy_reads: 4
distinct_valid_solutions_found: 1
multiplicity_evidence: unknown
shown: 1
solution 1:
534678912
672195348
198342567
859761423
426853791
713924856
961537284
287419635
345286179

500回の読み出しのうち、エネルギー0の状態は4回得られました。いずれも同じ完成盤面です。 まだ観測していない別解が存在する可能性を除けないため、この結果だけでは一意解とは判定できません。

解がない問題

fixtures/unsat-9x9.sdk では、同じ実行条件で得た最良エネルギーは4でした。 ゼロエネルギー状態は見つかっていません。ただし、有限回のサンプリングで見つからなかったことは 解なしの証明にならないため、結果は unknown です。

解が複数ある問題

fixtures/multiple-9x9.sdk では、固定後の変数は36個です。同じ実行条件で、 500回中498回がゼロエネルギーとなり、検証を通る異なる盤面が2個得られました。 二つの有効な盤面を得たので、複数解と確定できます。 ただし、観測した個数は全解数ではありません。

この方法で分かること

本定式化において、エネルギー0で共通検証器を通過する状態は数独の解に対応します。したがって、 該当状態が1個得られれば solved と出力します。異なる2個の状態が得られた場合は一意解でないことを 確定します。

一方、掲載したサンプラーで評価できない項目は以下の2点です。

• 一盤面のみが得られた場合の解の一意性

• ゼロエネルギーの盤面が見つからなかった場合の解なし判定

解なしを判定するには、たとえばQUBOの最適値が正であることを厳密に証明します。 最適値が0だと分かっても、一意性はまだ確定しません。一意性には、見つけた盤面以外の ゼロエネルギー解が存在しないことを確認する必要があります。 SATやSMTへ同じ数独を渡し、最初の盤面を除外して別解がないことを証明する方法も使えます。 これらの保証は、QUBOという定式化だけでなく、ソルバーと追加の検査によって決まります。

量子アニーラでも、ハードウェアから返った標本を評価し、最適性や解なしを確認する手順が必要です。 掲載結果はCPU上のnealで得たものです。

参考文献

[1]

Fred Glover, Gary Kochenberger, and Yu Du. Quantum Bridge Analytics I: A Tutorial on Formulating and Using QUBO Models. 4OR, 17(4):335–371, 2019. URL: https://arxiv.org/abs/1811.11538, doi:10.1007/s10288-019-00424-y.

[2]

D-Wave Systems Inc. dwave.samplers.SimulatedAnnealingSampler.sample. D-Wave Quantum Computing Products Documentation, 2026. Accessed 2026-08-08. URL: https://docs.dwavequantum.com/en/latest/ocean/api_ref_samplers/generated/dwave.samplers.SimulatedAnnealingSampler.sample.html.