第14章 反復射影で数独を解く¶
数独の候補を、0か1に決める前の実数として持ってみます。あるマスの四候補が
[0.1, 0.7, -0.2, 0.3] なら、値が0.7で最も大きい2番目の候補を1にし、残りを0にすると、
「一つのマスに数字が一つ」という規則を満たせます。今度は行の規則に合わせて候補を直すと、
先ほど決めたマスが崩れるかもしれません。それでも二種類の修正を繰り返し、すべての規則を
同時に満たす点を探します。
このように、点を制約集合内の最も近い点へ移す操作を繰り返す方法を 反復射影 と呼びます。
数独は1から9（または1から4）の整数を入れる離散問題ですが、反復射影では三つのステップで 扱います。まず各候補を0か1の二値ではなく実数へ拡張し、連続空間上の点として持ちます。次に 二種類の射影を繰り返し、すべての制約を同時に満たす点へ近づけます。最後に、各マスで最も 大きい候補を1、残りを0として離散的な盤面へ戻します。途中の実数値は、候補がどれだけ有力か を表す作業用の値であり、最終的な整数の答えそのものではありません。
二本の直線へ順番に射影する¶
最初に、平面上の点 \((x,y)\) を二本の直線へ移してみます。一つ目の集合 \(A\) は \(x\) 軸、二つ目の集合 \(B\) は直線 \(y=x\) です。共通部分は原点だけです。
点から集合内の最も近い点へ移す操作が 射影 です。\(x\) 軸への射影は \(y\) を0に します。直線 \(y=x\) への射影では、二つの座標を平均します。
def project_x_axis(point: Vector) -> Vector:
"""点をx軸上の最も近い位置へ移す。"""
return np.array([point[0], 0.0])
def project_diagonal(point: Vector) -> Vector:
"""点を直線y=x上の最も近い位置へ移す。"""
mean = (point[0] + point[1]) / 2.0
return np.array([mean, mean])\(P_A\) と \(P_B\) を順番に適用すると、点は次のように動きました。
$ uv run --with-requirements examples/13-iterative-projection/requirements.txt \
python examples/13-iterative-projection/small_projection.py
start: (2.000000, 1.000000)
iteration 1: (1.000000, 1.000000)
iteration 2: (0.500000, 0.500000)
iteration 3: (0.250000, 0.250000)
iteration 4: (0.125000, 0.125000)
iteration 5: (0.062500, 0.062500)
一回目の交互射影後は \((1,1)\)、二回目は \((0.5,0.5)\) です。原点に向かって座標が半分ずつに なっています。これは二本の直線についての結果です。一般の閉凸集合にも交互射影の収束理論が ありますが、数独で使う「候補の一つだけが1」という集合は凸集合ではありません。候補0と候補1を 結ぶ途中の点は、どちらの候補も0と1の間になり、数独の割り当てではないためです。非凸な集合では、 同じ二つの射影を使っても交点へ収束する保証はありません。[1]
局所制約を分けてから合意させる¶
数独には、マス、行、列、ブロックの規則があります。各規則への射影は個別に書けます。そこで、 すべてを一度に扱う代わりに「行 \(r\)、列 \(c\) に数字 \(d\) を置く候補」 \(x_{r,c,d}\) を四つ複製します。一つはマスの規則、残りは行、列、ブロックの規則に使います。
四種類の複製を別々に修正すれば、各局所規則を満たす点は構成されます。ただし、同じ候補についてマス側は1、 行側は0と主張することがあります。四つの値を平均して同じ値へそろえれば、今度は局所規則が 崩れます。局所規則を満たす集合を \(D\)、四つの複製が一致する集合を \(C\) とします。 それぞれへの射影を \(P_D\)、\(P_C\) と書きます。この分け方はDivide and Concurと 呼ばれます。[2] 詳細は 原論文の公開版 でも確認できます。
図の数値は一候補についての例です。\(P_C\) は四つを平均し、\(P_D\) は各複製を 担当する局所制約へ戻します。破線と実線は二種類の射影を区別しています。¶
4×4数独なら、元の候補配列は \(4\times4\times4=64\) 要素です。四種類に複製するので、 反復中は \(4\times4\times4\times4=256\) 個の実数を持ちます。配列の先頭の添字を、 マス、行、列、ブロックのどの規則を担当するかに使います。
軸 |
意味 |
4×4での値 |
|---|---|---|
|
担当する規則 |
マス、行、列、ブロックの4種類 |
|
行 |
0から3 |
|
列 |
0から3 |
|
数字の候補 |
0から3。盤面の数字1から4に対応 |
局所制約への射影¶
\(P_D\) は四つの複製を互いに相談させず、それぞれの局所制約だけを満たします。
- マスの複製
各マスで四候補の最大値を一つ選び、そこを1、残りを0にします。初期配置のあるマスは、 最大値に関係なく与えられた数字を選びます。
- 行の複製
行と数字を一組にし、その数字を置く列を一つ選びます。
- 列の複製
列と数字を一組にし、その数字を置く行を一つ選びます。
- ブロックの複製
2×2ブロックと数字を一組にし、その数字を置くマスを一つ選びます。
一要素だけが1で、残りが0の配列を one-hotベクトル と呼びます。四候補から最大値を選ぶ処理は、 実数の四次元ベクトルからone-hotベクトルへのユークリッド距離を 最小にする射影です。最大値が同じ候補が複数あれば最も近い点も複数あります。この実装では NumPyの argmax が最初に 見つけた候補を選びます。
def project_divide(values: FloatArray, givens: Board) -> FloatArray:
"""四種類の局所制約を、複製ごとに独立して満たす。"""
_, side, _, _ = values.shape
box_side = isqrt(side)
projected = np.zeros_like(values)
# 複製0では、各マスで最も大きな候補を一つだけ残す。
for row in range(side):
for column in range(side):
given = givens[row * side + column]
digit = given - 1 if given else int(np.argmax(values[0, row, column]))
projected[0, row, column, digit] = 1.0
# 複製1では、行ごとに各数字を置く列を一つだけ選ぶ。
for row in range(side):
for digit in range(side):
column = int(np.argmax(values[1, row, :, digit]))
projected[1, row, column, digit] = 1.0
# 複製2では、列ごとに各数字を置く行を一つだけ選ぶ。
for column in range(side):
for digit in range(side):
row = int(np.argmax(values[2, :, column, digit]))
projected[2, row, column, digit] = 1.0
# 複製3では、ブロックごとに各数字を置くマスを一つだけ選ぶ。
for top in range(0, side, box_side):
for left in range(0, side, box_side):
for digit in range(side):
block = values[
3,
top : top + box_side,
left : left + box_side,
digit,
]
offset = int(np.argmax(block))
row_offset, column_offset = divmod(offset, box_side)
projected[3, top + row_offset, left + column_offset, digit] = 1.0
return projected初期配置を固定するのはマスの複製だけです。残る三つの複製は、合意射影を通じて同じ値へ近づきます。 二つの集合の交点では四複製が一致するので、初期配置は行、列、ブロック側にも反映されています。
合意集合への射影¶
\(P_C\) は、同じ \(x_{r,c,d}\) に対応する四つの値を平均し、その平均を四複製すべてへ 書き戻します。
def project_concur(values: FloatArray) -> FloatArray:
"""同じ候補に対する四つの複製を平均し、同じ値にそろえる。"""
consensus = values.mean(axis=0, keepdims=True)
return np.broadcast_to(consensus, values.shape).copy()平均後の配列は、各候補について四つの値が同じなので \(C\) に属します。平均値が0か1である 必要はありません。合意することだけを受け持ち、数独の局所規則は \(P_D\) に任せています。
二つの射影をDifference Mapで組み合わせる¶
この実装では、単純に \(P_D\) と \(P_C\) を交互適用しません。現在の実数配列を \(x_k\) とし、まず \(D\) へ射影した点 \(d_k\) を求めます。\(d_k\) で \(x_k\) を反射した点 \(2d_k-x_k\) を \(C\) へ射影し、\(c_k\) とします。 反射とは、\(d_k\) を鏡の位置として、\(x_k\) を反対側へ同じ距離だけ移す操作です。
二つの射影結果の差を使い、次の配列を更新します。
これはDifference Mapの一形態で、relaxed reflect-reflect（RRR）とも呼ばれる更新です。 Difference Mapは、単純な交互射影で起きる停滞への対処として設計されていますが、非凸問題で 必ず交点を見つける保証はありません。[3]
停止判定には、\(c_k\) と \(d_k\) の各要素の差のうち、絶対値が最大のものを使います。 これをこの章の 残差 \(r_k\) とします。
残差が0なら \(c_k=d_k\) です。\(d_k\) は局所制約を満たす集合 \(D\) 上にあり、
\(c_k\) は四複製が合意する集合 \(C\) 上にあります。二点が一致すれば、その点は
\(D\cap C\) に属し、局所規則と合意を同時に満たします。そこから完成盤面を取り出せます。
計算機では丸め誤差があるため、
tolerance 以下で停止します。停止後にも盤面を共通検証器へ渡し、初期配置、行、列、ブロックを
すべて満たす場合だけ solved とします。
def solve(
givens: Board,
*,
seed: int,
max_iterations: int,
tolerance: float,
beta: float,
) -> Result:
side, _ = board_geometry(givens)
rng = np.random.default_rng(seed)
values = rng.normal(size=(4, side, side, side))
initial_residual = float("nan")
final_residual = float("inf")
for iteration in range(1, max_iterations + 1):
# 局所制約を満たす点で反射し、その点を合意集合へ射影する。
divided = project_divide(values, givens)
concurred = project_concur(2.0 * divided - values)
difference = concurred - divided
final_residual = float(np.max(np.abs(difference)))
if iteration == 1:
initial_residual = final_residual
# Difference Mapの一回分。残差0なら二つの射影結果が一致する。
values = values + beta * difference
if final_residual <= tolerance:
candidate = decode(divided)
valid, errors = validate_solution(candidate, givens)
if valid:
return Result(
status="solved",
seed=seed,
beta=beta,
tolerance=tolerance,
max_iterations=max_iterations,
iterations=iteration,
initial_residual=initial_residual,
final_residual=final_residual,
solution=candidate,
)
return Result(
status="unknown",
seed=seed,
beta=beta,
tolerance=tolerance,
max_iterations=max_iterations,
iterations=iteration,
initial_residual=initial_residual,
final_residual=final_residual,
solution=None,
validation_errors=errors,
)
# 反復上限までに交点へ到達しなくても、解なしとは結論しない。
return Result(
status="unknown",
seed=seed,
beta=beta,
tolerance=tolerance,
max_iterations=max_iterations,
iterations=max_iterations,
initial_residual=initial_residual,
final_residual=final_residual,
solution=None,
)beta は一回の更新量を調整します。掲載結果では0.5に固定しました。初期配列は標準正規分布から
作るため、seed も結果の一部です。反復上限到達時および残差が閾値以下であっても復元盤面が
検証失敗となった不整合状態においては、一律 unknown として処理します。
4×4問題を固定条件で実行する¶
詳しい実行手順やオプションについては examples/13-iterative-projection/README.md を参照してください。
固定した初期シードとパラメータ（\(\beta=0.5\) 等）での試行において、43反復で残差が許容誤差（\(10^{-8}\) 以下）まで減少し、検証を通過する有効な完成盤面が得られます。ただし、この反復数は初期値やパラメータ設定に依存した特定試行の結果です。
収束しなければunknownにする¶
解が存在しない盤面（あるいは特定の初期値において収束しなかった盤面）に対して実行すると、最大反復上限（2,000反復など）に達しても残差が許容誤差以下まで低下しません。
この場合、結果は unsat ではなく unknown となります。反復射影における交点未検出は制約集合の交点非存在（解なし）を厳密に証明するものではなく、単に指定条件で収束しなかった状態を表すためです。
また、複数解を持つ問題で一解が得られた場合であっても、Difference Mapの1回の試行のみから他の解の存在や一意性を確定することはできません。連続空間での反復探索は解の列挙を行わないためです。
この方法で分かること¶
solved の盤面は、残差による停止判定に加えて共通検証器を通しています。そのため、返した盤面が
数独の規則と初期配置を満たすことは確認できます。途中の実数配列や、残差が小さくなった事実
だけを解として扱ってはいません。
unknown は、指定したシード、パラメータ、反復上限で検証済み盤面を得られなかった状態を示します。
この結果から解の有無や一意性は判定しません。初期配列の生成以外で乱数は使用せず、更新プロセスは
決定論的です。シードを変更した場合は反復数、解出力、収束挙動も変化します。
反復射影では、離散的な数独を実数空間の二つの集合の交点として扱います。射影を個別に記述できれば、 同じ枠組みをほかの制約充足問題にも適用します。[2] 掲載実装は射影で 見つけた候補を検証して採用するヒューリスティックであり、解なしや一意性の判定には、バックトラック、 SAT、SMT、BDD、整数計画法などの厳密解法が必要です。
参考文献¶
Eric C. Chi and Kenneth Lange. Techniques for Solving Sudoku Puzzles. 2012. URL: https://arxiv.org/abs/1203.2295, arXiv:1203.2295, doi:10.48550/arXiv.1203.2295.
Simon Gravel and Veit Elser. Divide and Concur: A General Approach to Constraint Satisfaction. Physical Review E, 78(3):036706, 2008. URL: https://arxiv.org/abs/0801.0222, doi:10.1103/PhysRevE.78.036706.
Veit Elser. Phase Retrieval by Iterated Projections. Journal of the Optical Society of America A, 20(1):40–55, 2003. URL: https://opg.optica.org/josaa/abstract.cfm?uri=josaa-20-1-40, doi:10.1364/JOSAA.20.000040.