第14章反復射影で数独を解く¶

数独の候補を、0か1に決める前の実数として持ってみます。あるマスの九候補が [0.1, 0.7, -0.2, 0.3, 0.0, -0.1, 0.2, 0.4, 0.1] なら、値が0.7で最も大きい2番目の候補を1にし、残りを0にすると、「一つのマスに数字が一つ」という規則を満たせます。今度は行の規則に合わせて候補を直すと、先ほど決めたマスが崩れるかもしれません。それでも二種類の修正を繰り返し、すべての規則を同時に満たす点を探します。

このように、点を制約集合内の最も近い点へ移す操作を繰り返す方法を反復射影と呼びます。

数独は各マスに離散的な整数を割り当てる問題ですが、反復射影では候補を0と1の間に限らない実数として扱い、連続空間上の幾何学的な点として探索します。マスや行などの局所的な制約を満たす点と、複製同士が一致する合意集合という二つの集合を行き来しながら、交点へ近づけていきます。最後に各マスで最も値が大きい候補を選んで整数の盤面へ復元します。反復中の実数値は交点を探すための作業用の座標であり、完成盤面の数値そのものとは異なります。

二本の直線へ順番に射影する¶

最初に、平面上の点 \((x,y)\) を二本の直線へ移してみます。一つ目の集合 \(A\) は \(x\) 軸、二つ目の集合 \(B\) は直線 \(y=x\) です。共通部分は原点だけです。

点から集合内の最も近い点へ移す操作が射影です。\(x\) 軸への射影は \(y\) を0にします。直線 \(y=x\) への射影では、二つの座標を平均します。

\[P_A(x,y)=(x,0), \qquad P_B(x,y)=\left(\frac{x+y}{2},\frac{x+y}{2}\right)\]

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\) は各複製を担当する局所制約へ戻します。破線と実線は二種類の射影を区別しています。¶

9×9数独なら、元の候補配列は \(9\times9\times9=729\) 要素です。四種類に複製するので、反復中は \(4\times9\times9\times9=2{,}916\) 個の実数を持ちます。配列の先頭の添字を、マス、行、列、ブロックのどの規則を担当するかに使います。

候補配列 `values[k, r, c, d]` の四つの軸¶
軸	意味	9×9での値
`k`	担当する規則	マス、行、列、ブロックの4種類
`r`	行	0から8
`c`	列	0から8
`d`	数字の候補	0から8。盤面の数字1から9に対応

局所制約への射影¶

\(P_D\) は四つの複製を互いに相談させず、それぞれの局所制約だけを満たします。

マスの複製: 各マスで九候補の最大値を一つ選び、そこを1、残りを0にします。初期配置のあるマスは、最大値に関係なく与えられた数字を選びます。
行の複製: 行と数字を一組にし、その数字を置く列を一つ選びます。
列の複製: 列と数字を一組にし、その数字を置く行を一つ選びます。
ブロックの複製: 3×3ブロックと数字を一組にし、その数字を置くマスを一つ選びます。

一要素だけが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}\) に対応する四つの値を平均し、その平均を四複製すべてへ書き戻します。

\[\bar{x}_{r,c,d}=\frac{1}{4} \sum_{k=0}^{3}x_{k,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\) を反対側へ同じ距離だけ移す操作です。

\[d_k=P_D(x_k), \qquad c_k=P_C(2d_k-x_k)\]

二つの射影結果の差を使い、次の配列を更新します。

\[x_{k+1}=x_k+\beta(c_k-d_k)\]

これはDifference Mapの一形態で、relaxed reflect-reflect（RRR）とも呼ばれる更新です。 Difference Mapは、単純な交互射影で起きる停滞への対処として設計されていますが、非凸問題で必ず交点を見つける保証はありません。[3]

停止判定には、\(c_k\) と \(d_k\) の各要素の差のうち、絶対値が最大のものを使います。これをこの章の残差 \(r_k\) とします。

\[r_k=\lVert c_k-d_k\rVert_\infty\]

残差が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 とします。残差が閾値以下でも、復元した盤面が検証を通らなければ unknown です。

9×9問題を固定条件で実行する¶

詳しい実行手順やオプションについては examples/13-iterative-projection/README.md を参照してください。

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

$ uv run --frozen python examples/13-iterative-projection/solve.py \
    fixtures/standard-9x9.sdk --seed 0 --max-iterations 2000 \
    --beta 0.5 --tolerance 1e-8

status: solved
seed: 0
beta: 0.5
tolerance: 1e-08
max_iterations: 2000
iterations: 1528
initial_residual: 1.69227635
final_residual: 5.83694926e-09
solution:
534678912
672195348
198342567
859761423
426853791
713924856
961537284
287419635
345286179

この条件では1,528反復で残差が \(10^{-8}\) 以下になり、共通検証器を通る盤面が得られました。 1,528回という反復数は、この入力、シード、パラメータでの結果です。

収束しなければunknownにする¶

入力を fixtures/unsat-9x9.sdk に替え、同じパラメータで実行すると、 2,000反復しても残差が許容値まで下がりませんでした。

status: unknown
seed: 0
beta: 0.5
tolerance: 1e-08
max_iterations: 2000
iterations: 2000
initial_residual: 1.69227635
final_residual: 1.16552723

結果は unknown です。指定条件で交点を見つけられなかっただけでは、交点が存在しないとは証明できません。解がある盤面でも、初期値や反復上限によっては unknown になります。

また、複数解を持つ問題で一解が得られた場合であっても、Difference Mapの1回の試行のみから他の解の存在や一意性を確定することはできません。連続空間での反復探索は解の列挙を行わないためです。

この方法で分かること¶

solved の盤面は、残差による停止判定に加えて共通検証器を通しています。そのため、返した盤面が数独の規則と初期配置を満たすことは確認できます。途中の実数配列や、残差が小さくなった事実だけを解として扱ってはいません。

unknown は、指定したシード、パラメータ、反復上限で検証済み盤面を得られなかった状態を示します。この結果から解の有無や一意性は判定しません。初期配列の生成以外で乱数は使用せず、更新プロセスは決定論的です。シードを変更した場合は反復数、解出力、収束挙動も変化します。

反復射影では、離散的な数独を実数空間の二つの集合の交点として扱います。射影を個別に記述できれば、同じ枠組みをほかの制約充足問題にも適用します。[2] 掲載実装は射影で見つけた候補を検証して採用するヒューリスティックであり、解なしや一意性の判定には、バックトラック、 SAT、SMT、BDD、整数計画法などの厳密解法が必要です。

参考文献¶

[1]

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.

[2] (1,2)

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.

[3]

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.

第14章 反復射影で数独を解く¶