第15章 因子グラフで数独を解く¶
空いている一マスについて、数字1から4の重みを [0.1, 0.6, 0.2, 0.1] のように持つことに
します。このマスは同じ行、列、ブロックにいるマスと関係しています。行の側から「2はほかの
マスと両立しにくい」、列の側から「3なら両立しやすい」という情報を受け取って重みを更新します。
今度は更新後の重みを周囲へ返します。このやり取りを繰り返し、各マスで最も重い数字を選びます。
ここではまず、重みを候補同士の有力さを比べる非負の数として読みます。ループのある数独では、
正規化して合計を1にしても、厳密な周辺確率になるとは限りません。
数独は1から4の整数を入れる離散問題ですが、因子グラフでは三つのステップで扱います。まず 各候補に非負の実数の重みを割り当て、数字の有力さを連続的な値として持ちます。次に変数と 因子の間で重みを繰り返し交換し、周囲の情報を反映させます。最後に、各マスで最も重い候補の 数字を選んで離散的な盤面へ戻します。途中の重みは候補間の相対的な有力さを表す近似値であり、 最終的な整数の答えそのものではありません。
この関係を図にしたものが 因子グラフ です。マスを表す 変数ノード と、 行、列、ブロックの規則を表す 因子ノード を別の形で描きます。辺で直接つながるのは、変数と その変数を使う因子だけです。ノード間で重みの配列を交換する計算を belief propagation と 呼びます。因子グラフとメッセージ計算は、誤り訂正符号の復号や確率的な推論にも使われています。 [1]
変数と規則をノードにする¶
16マスを \(X_1,\ldots,X_{16}\) とし、各変数は1から4を 取ります。行、列、2×2ブロックはそれぞれ四つずつあるので、因子ノードは12個です。一つのマスは 行、列、ブロックの三因子につながり、辺は全部で48本になります。
因子は、つながっている四変数がすべて異なれば1、同じ数字があれば0を返します。この関数を \(\psi_u\) とします。初期配置を表す関数 \(\phi_i\) は、指定された数字だけを1とし、 それ以外を0にします。空きマスではどの数字も1です。盤面全体の関数は、これらの積で書けます。
\(F(\boldsymbol{x})=1\) になる割り当てが数独の解です。大きな関数を、少数の変数だけを見る 小さな関数の積に分けた形が、因子グラフの辺に対応します。
円は変数ノード、四角は因子ノードです。一つの辺では両方向に別のメッセージを送ります。 数独では行因子と列因子だけを取り出しても閉じた経路ができます。¶
メッセージが運ぶもの¶
一つのメッセージは、数字ごとの非負の重みです。変数 \(X\) から因子 \(f\) への メッセージでは、初期配置の重みと、\(f\) 以外の隣接因子から届いた重みを掛けます。 次の式の \(\propto\) は、右辺を計算した後、全候補の合計が1になるよう同じ倍率で割ることを 表します。
\(N(X)\) は \(X\) につながる因子の集合です。送り先 \(f\) から受け取った値を掛け ないのは、同じ情報をすぐ送り返して二重に数えないためです。
因子から変数へのメッセージでは、周囲の変数を規則どおりに割り当てた場合を調べます。この実装では、 それらの重みを合計する sum-product と、最大の重みだけを選ぶ max-product の二方式を使います。 この因子では、1から4の順列が「すべて異なる」割り当てです。sum-productでは、対象のマスを数字 \(d\) にした各割り当ての重みを足します。
次の式の \(f(\boldsymbol{x})\) は、因子ノード \(f\) が表す関数で、前節の \(\psi_u\) に当たります。\(N(f)\) は、\(f\) につながる変数の集合です。
max-productでは式の和を最大値に置き換えます。どちらも計算後に合計が1になるよう正規化します。
一つの因子だけで計算する¶
まず、三変数が1から3を一度ずつ取る因子を試します。対象を \(X\) とし、残る二変数から
[0.6, 0.3, 0.1] と [0.2, 0.3, 0.5] が届いたとします。コードは六通りの順列を列挙し、
対象以外の重みを掛けます。
def all_different_message(
incoming: list[Message] | tuple[Message, ...],
target_position: int,
method: Method,
) -> Message | None:
"""all-different因子から一つの変数へ送るメッセージを計算します。"""
side = len(incoming)
outgoing = [0.0] * side
# 一つの順列が、因子を満たす数字の割り当て一通りに対応します。
for assignment in permutations(range(side)):
weight = prod(
incoming[position][assignment[position]]
for position in range(side)
if position != target_position
)
target_digit = assignment[target_position]
if method == "sum-product":
# 周囲の全割り当てから届く重みを足し合わせます。
outgoing[target_digit] += weight
else:
# 周囲の割り当てのうち、最大の重みだけを残します。
outgoing[target_digit] = max(outgoing[target_digit], weight)
return normalize(outgoing)#!/usr/bin/env python3
"""Compare sum-product and max-product messages for one factor."""
from __future__ import annotations
from solve import all_different_message
def show(values: tuple[float, ...] | None) -> str:
if values is None:
return "none"
return "[" + ", ".join(f"{value:.3f}" for value in values) + "]"
def main() -> None:
# xへのメッセージを求めるので、x自身から届いた値は計算に使われません。
incoming = (
(1 / 3, 1 / 3, 1 / 3),
(0.6, 0.3, 0.1),
(0.2, 0.3, 0.5),
)
print("sum-product:", show(all_different_message(incoming, 0, "sum-product")))
print("max-product:", show(all_different_message(incoming, 0, "max-product")))
if __name__ == "__main__":
main()
sum-product: [0.243, 0.432, 0.324]
max-product: [0.238, 0.476, 0.286]
sum-productで \(X=2\) とする場合を手で追うと、残る二変数には1と3を異なる順で割り当てます。
正規化前の重みは \(0.6\times0.5+0.1\times0.2=0.32\) です。同様に \(X=1,3\) の重みは
0.18と0.24なので、合計0.74で割ると [0.243, 0.432, 0.324] になります。コードの
permutations は割り当てを列挙し、prod は各割り当ての重みを掛けています。
二つの方式とも数字2を最も重くしていますが、三つの重みは一致しません。sum-productでは有効な 割り当てを足し、max-productでは最大の割り当てだけを使うためです。この差は、因子グラフの形では なく、メッセージをまとめる演算から生じます。
数独のグラフで反復する¶
実装は、初期配置のマスを指定数字だけが1となるone-hotの重み、空きマスを一様な重みで初期化します。 因子から変数へ送る重みも最初は一様です。一回の反復で、まず全因子から全変数へのメッセージを更新し、 その値を使って全変数から全因子へのメッセージを更新します。
数独のグラフには短いループがあります。同じ情報が何周もして重みが振動するのを抑えるため、更新値
だけへ一度に置き換えず、前回の値を damping=0.5 の割合で残します。全メッセージについて前回
からの最大変化を 残差 とし、残差が \(10^{-10}\) 未満なら停止します。反復上限は200回です。
residual = float("inf")
for iteration in range(1, max_iterations + 1):
new_factor_to_variable = [[uniform for _ in unit] for unit in factors]
residual = 0.0
# 各因子は、同じ単位にいるほかのマスの重みをまとめます。
for factor_id, unit in enumerate(factors):
incoming = variable_to_factor[factor_id]
for position, variable in enumerate(unit):
proposed = all_different_message(incoming, position, method)
if proposed is None:
return SolveResult(
"unknown",
method,
iteration,
False,
residual,
0,
reason_unknown="因子から送れる正の重みがなくなりました",
)
old = factor_to_variable[factor_id][position]
# ループで値が振動しにくいよう、前回の値を半分残します。
message = damp(old, proposed, damping)
new_factor_to_variable[factor_id][position] = message
residual = max(residual, difference(old, message))
new_variable_to_factor = [[uniform for _ in unit] for unit in factors]
# 各マスは、送り先以外の因子から届いた重みを掛け合わせます。
for variable, edges in enumerate(neighbors):
for target_factor, target_position in edges:
weights = [
priors[variable][digit]
* prod(
new_factor_to_variable[factor_id][position][digit]
for factor_id, position in edges
if factor_id != target_factor
)
for digit in range(side)
]
proposed = normalize(weights)
if proposed is None:
return SolveResult(
"unknown",
method,
iteration,
False,
residual,
0,
reason_unknown="変数から送れる正の重みがなくなりました",
)
old = variable_to_factor[target_factor][target_position]
message = damp(old, proposed, damping)
new_variable_to_factor[target_factor][target_position] = message
residual = max(residual, difference(old, message))
factor_to_variable = new_factor_to_variable
variable_to_factor = new_variable_to_factor
if residual < tolerance:
break
else:
return SolveResult(
"unknown",
method,
max_iterations,
False,
residual,
0,
reason_unknown="反復上限までにメッセージが収束しませんでした",
)
最後に、各変数について初期値と三因子から届いた重みを掛け、正規化します。この数字ごとの重みを
belief と呼びます。最大値となる候補が複数あって同点の場合、あるいは上位2つの差が \(10^{-7}\) 以内と僅差の場合には、数字を確定させず unknown
にします。一意な最大値を選べても、独立に選んだ16個の数字が数独の規則を満たすとは限りません。
完成盤面を共通検証器へ渡し、初期配置、行、列、ブロックをすべて確認できた場合だけ solved
とします。別の探索や盤面の修復はしません。
# 各変数のbeliefは、初期値と隣接因子から届いた全メッセージの積です。
beliefs: list[Message] = []
for variable, edges in enumerate(neighbors):
weights = [
priors[variable][digit]
* prod(
factor_to_variable[factor_id][position][digit]
for factor_id, position in edges
)
for digit in range(side)
]
belief = normalize(weights)
if belief is None:
return SolveResult(
"unknown",
method,
iteration,
True,
residual,
0,
reason_unknown="beliefを正規化できませんでした",
)
beliefs.append(belief)
solution, ties, reason = decode(tuple(beliefs), board, tie_tolerance)
if solution is None:
return SolveResult(
"unknown",
method,
iteration,
True,
residual,
ties,
reason_unknown=reason,
)
return SolveResult("solved", method, iteration, True, residual, ties, solution)
木では厳密、数独では近似になる¶
因子グラフにループがなければ、端から内側へメッセージを送り、逆向きにも送ることで、sum-productは 各変数の厳密な周辺値を計算できます。max-productも同じ木構造上で最大重みの割り当てを厳密に 求められます。[1]
数独の因子グラフには、図の右側のようなループが多数あります。そこでも同じ更新式を反復する方法を loopy belief propagationと呼びます。メッセージが収束する保証はなく、収束しても厳密な周辺値とは 限りません。[2] 数独へ適用した研究でも、短いループを含むグラフで 反復し、すべての問題を解ける方式ではないことが報告されています。[3]
この章の belief は、ループ付きグラフ上で数字を選ぶための近似的な重みとして扱います。残差が
小さくなったことは、メッセージが固定点に近づいたという意味です。それだけで、復元した盤面が正しい
とも、ほかの解がないとも言えません。そのため停止判定とは別に盤面を検証します。
4×4で実装する理由¶
このサンプルでは、因子から一マスへメッセージを送るたびに数字の順列をすべて列挙します。4×4なら 一因子につき \(4!=24\) 通りですが、9×9では \(9!=362{,}880\) 通りです。それを各因子の 各辺で毎反復するため、9×9の教材用実装としては無駄が大きくなります。
4×4を使うのは因子グラフ自体の制限ではありません。順列の全列挙を避ける動的計画法や、制約の分け方 を変えた実装も考えられます。この章では、メッセージの意味をコードから直接読めることを優先し、 4×4だけを明示的に受け付けます。
実行と検証¶
実行手順とオプションは examples/14-factor-graph/README.md にまとめています。同じ一意解問題を sum-productとmax-productで解き、複数解問題と矛盾問題をsum-productで調べた結果は次のとおりです。
問題 |
方式 |
|
反復数 |
残差 |
復元結果 |
|---|---|---|---|---|---|
一意解 |
sum-product |
|
44 |
|
検証済み盤面、同率0マス |
一意解 |
max-product |
|
43 |
|
検証済み盤面、同率0マス |
複数解 |
sum-product |
|
44 |
|
最大beliefが同率のマス8個 |
矛盾 |
sum-product |
|
200 |
|
反復上限で未収束 |
一意解問題では両方式とも同じ4×4盤面を返し、共通検証器を通りました。sum-productは44回、
max-productは43回で収束しています。ただし、一盤面を得ただけではほかの解を除外していないため、
どちらも unique: undetermined です。
複数解問題では数字を決められない¶
複数解を持つ問題に対してbelief propagationを実行すると、メッセージは収束するものの、複数のマスで候補の近似的な重み（belief）が同率となる拮抗状態が発生します。バックトラック等の分岐復元処理を含まないため、この場合は解を確定できず unknown となります。
矛盾問題もunknownにする¶
矛盾問題（解なし）に対して実行した場合、ループ付きグラフ上のメッセージ更新が収束せず、最大反復上限に達して結果は unknown となります。loopy belief propagation におけるメッセージの未収束は、解の非存在を証明する根拠にはならないためです。
この方法で分かること¶
この実装が solved と返す盤面は、共通検証器が数独の規則と初期配置を確認しています。近似的な
重みから作った盤面でも、検証を通った完成盤面そのものは有効な解です。
unknown には、メッセージが収束しなかった場合、最大候補が同率の場合、復元盤面が規則を破った
場合が含まれます。どの場合も解なしを意味しません。一意性を調べる全解列挙や、最初の解を除外して
再探索する機構もないため、結果は常に unique: undetermined です。
局所的な規則同士で重みを交換し、検証を通る盤面へ到達する場合があります。ループ付きグラフの収束は、 数独の解の存在や一意性を証明しません。そこまで判定するには、バックトラック、SAT、SMT、BDD、 整数計画法などの厳密解法を用います。
参考文献¶
Frank R. Kschischang, Brendan J. Frey, and Hans-Andrea Loeliger. Factor Graphs and the Sum-Product Algorithm. IEEE Transactions on Information Theory, 47(2):498–519, 2001. URL: https://www.isiweb.ee.ethz.ch/papers/arch/aloe-2001-1.pdf, doi:10.1109/18.910572.
Alexander T. Ihler, John W. Fisher, III, and Alan S. Willsky. Loopy Belief Propagation: Convergence and Effects of Message Errors. Journal of Machine Learning Research, 6(31):905–936, 2005. URL: https://jmlr.org/papers/v6/ihler05a.html.
Todd K. Moon and Jacob H. Gunther. Multiple Constraint Satisfaction by Belief Propagation: An Example Using Sudoku. In 2006 IEEE Mountain Workshop on Adaptive and Learning Systems, 122–126. 2006. doi:10.1109/SMCALS.2006.250702.