第13章 QUBOで数独を解く¶
数独の規則を、守るか破るかで点数が変わる式にしてみます。正しい完成盤面なら0点、同じ行に 同じ数字を置くなどの違反があれば加点される式です。あとは、この点数を最小にする0と1の組合せを 探します。
このようなモデルを QUBO（Quadratic Unconstrained Binary Optimization） と呼びます。日本語にすれば「制約なし二値二次最適化」です。QUBOは問題の書き方であり、 解法や計算機を指定しません。掲載コードでは、通常のCPU上で焼きなまし法を実行します。
一つだけ選ぶ条件を点数にする¶
まず、二つの候補 \(x_1,x_2\) から一つだけ選ぶ場面を考えます。選ぶ候補を1、選ばない候補を 0とすれば、次の式を罰点にできます。
一つだけ選んだときは括弧の中が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になります。
一般に、候補 \(x_1,\ldots,x_n\) から一つだけ選ぶ条件も、同じ形で書けます。
右辺には、各変数だけに掛かる一次の項と、二変数の組に掛かる二次の項しかありません。条件を
ソルバーへ別に渡す代わりに、破ったときに大きくなる式へ埋め込んでいます。このため
unconstrained と呼ばれます。ただし、変数が0か1である制限は残っています。
QUBOは何を最小化するか¶
QUBOでは、二値変数のベクトル \(\boldsymbol{x}\) に対し、一次と二次の項を足したエネルギー \(E(\boldsymbol{x})\) を最小化します。ここでいうエネルギーは、候補の良し悪しを比べる点数です。 低いほど罰点が少なくなります。この章では係数の意味が見えるよう、次の書き方をします。
\(C\) は定数、\(a_i\) は一変数の係数、\(b_{ij}\) は二変数の組の係数です。文献や
ライブラリでは行列 \(Q\) を使って \(\boldsymbol{x}^{\mathsf T}Q\boldsymbol{x}\) と書く
こともあります。対称な成分をどう数えるかで係数の置き方が変わるため、コードでは
dimod.BinaryQuadraticModel の一次係数、二次係数、定数項を分けて渡します。
QUBOは、割当て、巡回路、スケジューリングなどの組合せ最適化を二値変数と罰則へ変換するときにも 使われます。制約を二次式へ移す手順と多くの定式化例は、Gloverらのチュートリアルにまとめられて います [1]。
QUBOは、古典的な厳密ソルバー、ヒューリスティック、量子アニーリングのいずれでも扱えます。 モデルを作った後、最低エネルギーを見つける処理は選んだソルバーが担当します。
数独を64個の二値変数で表す¶
この章では、9×9ではなく4×4の数独を使います。各ブロックは2×2で、各行、列、ブロックに1から4が 一度ずつ入ります。行 \(r\)、列 \(c\) のマスに数字 \(d\) を入れる候補を \(x_{r,c,d}\) とします。
4×4なら候補変数は \(4\times4\times4=64\) 個です。9×9では729個になり、変数間の二次項も 大幅に増えます。今回はサンプラーを短時間で何度も動かし、複数の低エネルギー状態を観察することを 優先して4×4にしました。盤面の大きさは違っても、候補変数と次の四種類の条件は同じです。
一つのマスには、数字が一つだけ入ります。
一つの行には、各数字が一度だけ現れます。
一つの列には、各数字が一度だけ現れます。
一つの2×2ブロックには、各数字が一度だけ現れます。
各条件は「四候補から一つだけ選ぶ」ので、先ほどの二次罰則をそのまま使えます。コードでは
add_exactly_one が、一次係数へ \(-1\)、候補ペアの二次係数へ \(2\)、定数項へ
\(1\) を加えます。weight は罰則全体に掛ける正の重みです。
たとえば左上のマスには、
\((x_{1,1,1}+x_{1,1,2}+x_{1,1,3}+x_{1,1,4}-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("この例は4×4数独だけを扱います")
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,
)
# 各2×2ブロックでも、それぞれの数字を一度だけ選ぶ。
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では、制約を破って目的値だけを小さくする状態に負けない よう、罰則の重みを決める必要があります。数独ではその競合がないため、重み選びの問題を避けられます。
焼きなましで状態をサンプリングする¶
モデルの構築には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)
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
board = decode_sample(datum.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)),
)焼きなまし法は、限られた反復回数で必ず最小エネルギーへ到達することを保証しません。サンプリング実行において
エネルギー0の状態に到達しなかった場合、プログラムは unsat ではなく unknown を返します。
「解の非存在」と「反復探索における未発見状態」を識別するためです。固定シードの設定は
掲載結果の再現性を確保しますが、最適性の厳密な証明を提供するものではありません。
実行と検証¶
モデルの実行および検証手順については examples/12-qubo/README.md を参照してください。
通常問題¶
サンプリング実行においてエネルギー0（全制約クリア）の状態が複数回観測されたものの、検証により得られた相異なる完成盤面が一つのみであるケースでは、確率的な重複観測にとどまり、この出力のみから一意解を厳密に確定することはできません（未観測解の存在を否定できないため）。
解が複数ある問題¶
空の4×4盤面などのサンプリングでは、エネルギー0となる異なる完成盤面が複数個得られます。互いに異なる二つ以上の有効な盤面が検証を通った時点で「複数解が存在する」と確定できます。ただし、観測された解の個数を全解数とみなすことはできません。
解がない問題¶
解が存在しない盤面では、アニーリングを実行しても最良エネルギーが1以上にとどまり、ゼロエネルギー状態を得られません。しかし、有限回のヒューリスティックサンプリングのみでは解の非存在を証明できないため、結果は unknown となります。
この方法で分かること¶
本定式化において、エネルギー0で共通検証器を通過する状態は数独の解に対応します。したがって、
該当状態が1個得られれば solved と出力します。異なる2個の状態が得られた場合は一意解でないことを
確定します。
一方、掲載したサンプラーで評価できない項目は以下の2点です。
一盤面のみが得られた場合の解の一意性
ゼロエネルギー状態未検出時における問題の不成立性（unsat）
これらを判定するには、QUBOを厳密に最適化して最適値を証明する手法や、前章までに扱ったSATや SMTのような完全な解法が必要です。掲載実装での限界は近似的な焼きなましサンプラーの採用による ものであり、同一のQUBOを完全解法ソルバーへ渡した場合は判定結果の性質が変化します。
量子アニーラでも、ハードウェアから返った標本を評価し、最適性や解なしを確認する手順が必要です。 掲載結果はCPU上のnealで得たものです。
参考文献¶
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.
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.