グラフの同型性判定という問題がある。ノードの数、エッジの数が同じ2つのグラフがあるときに、この2つが同じ連結構造をしているか、という問題。例えば、
A-B-C-D
と
P-Q-R
|
S
は、どのようなマッピングを行っても対応がつかないので、同型ではない。
整数計画法のソルバを使って解きたかったのだけど、うまく表現する方法が思いつかず。色々調べてみたらこちらを見つけた。
量子アニーラを前提としているようで、モデルがやや違うようだけど表現方法は同じでよさそうだったのでこれで書いてみた。
基本的な方針
基本的には、2つのグラフのノード列を縦横に並べた0/1の行列で表現する。たとえば以下はAがPにBがQにくるマッピングを表す。
| P | Q
--+---+---
A | 1 | 0
B | 0 | 1
この行列の各行、各列の和は1になる。これは、あるノードがマッピングする相手は1つしかない、ということを考えれば当然だ。
この条件だけではもちろんだめで、エッジが対応していなければならない。条件としては、
- グラフ1のノードu, vとグラフ2のノードs, t がそれぞれ対応しているなら、
- u, v間にエッジがあれば、s, t間にもエッジがなければならず
- u, v間にエッジがないなら、s, t間にもエッジはなってはならない
これらの条件をそのまま書こうとすると非常に大変。というか、書き方がわからない。ので、これを書き換える
- u, v 間にエッジがあり、s, t間にエッジが無いならば、u,v と s,t が対応していてはならない
- u, v 間にエッジがなく、s, t間にエッジがあるならば、u,v と s,t が対応していてはならない
これらが本当に等価の条件になっているのかよくわからないが、この方針で書いてみる。
u,v と s,t
コード
pyomoで書く。だんだん書き慣れてきたぞ。
from pyomo.environ import *
def get_v(edges):
s = set()
for i,j in edges:
s.add(i)
s.add(j)
return sorted(list(s))
def cross(X,Y):
return [(i,j) for i in X for j in Y]
def sort_order(tuples):
return [(i,j) if i <= j else (j, i) for i,j in tuples]
X = sort_order([("a", "b"), ("a", "c"), ("c", "d")])
Y = sort_order([("p", "q"), ("p", "r"), ("r", "s")])
X_v = get_v(X)
Y_v = get_v(Y)
model = ConcreteModel()
model.X_v = Set(initialize=X_v)
model.Y_v = Set(initialize=Y_v)
model.X_v.pprint()
model.Y_v.pprint()
# 割当のインデックス
model.map = Set(within=model.X_v * model.Y_v, initialize=cross(X_v, Y_v))
model.map.pprint()
# 割当変数 交差しているところが割当を表す
model.q = Var(model.map, within=Binary)
# 割当がユニークである制約
model.q_x = Constraint(model.Y_v, rule=lambda m, j: sum(m.q[(i,j)] for i in X_v) == 1)
model.q_y = Constraint(model.X_v, rule=lambda m, i: sum(m.q[(i,j)] for j in Y_v) == 1)
# Xにエッジu,vがあり場合 Yにエッジs,tがなければ、q[u,s]q[v,t]は 0
def exist_rule(m, u, v, s, t):
if u == v or s >= t:
return Constraint.Skip
if ((u, v) in X or (v, u) in X) and (s, t) not in Y:
return m.q[(u, s)] + m.q[(v, t)] <= 1
return Constraint.Skip
model.exist_edge = Constraint(model.X_v, model.X_v, model.Y_v, model.Y_v, rule=exist_rule)
# Yにエッジs,tがあり場合 Xにエッジu,vがなければ、q[u,s]q[v,t]は 0
def rexist_rule(m, u, v, s, t):
if u == v or s >= t:
return Constraint.Skip
if ((u, v) not in X or (u, v) not in X) and (s, t) in Y:
return m.q[(u, s)] + m.q[(v, t)] <= 1
return Constraint.Skip
model.rexist_edge = Constraint(model.X_v, model.X_v, model.Y_v, model.Y_v, rule=rexist_rule)
opt = SolverFactory('glpk')
opt.solve(model, tee=False)
model.display()
model.q.pprint()
ポイントはexist_ruleとrexist_rule だがこれであってるのか確信はない。
return m.q[(u, s)] + m.q[(v, t)] <= 1
の部分は、 uとs、vとtが同時に対応していることはない(どちらかが対応していてもOK)という条件。もとの記事だと
return m.q[(u, s)] * m.q[(v, t)] == 0
のように積で書かれていたが、線形計画法ではこのような式は書けなかったので上のように書いた。正しいんだろうか。
結果
上のコードを実行すると以下のように出力される。
Variables:
q : Size=16, Index=map
Key : Lower : Value : Upper : Fixed : Stale : Domain
('a', 'p') : 0 : 1.0 : 1 : False : False : Binary
('a', 'q') : 0 : 0.0 : 1 : False : False : Binary
('a', 'r') : 0 : 0.0 : 1 : False : False : Binary
('a', 's') : 0 : 0.0 : 1 : False : False : Binary
('b', 'p') : 0 : 0.0 : 1 : False : False : Binary
('b', 'q') : 0 : 1.0 : 1 : False : False : Binary
('b', 'r') : 0 : 0.0 : 1 : False : False : Binary
('b', 's') : 0 : 0.0 : 1 : False : False : Binary
('c', 'p') : 0 : 0.0 : 1 : False : False : Binary
('c', 'q') : 0 : 0.0 : 1 : False : False : Binary
('c', 'r') : 0 : 1.0 : 1 : False : False : Binary
('c', 's') : 0 : 0.0 : 1 : False : False : Binary
('d', 'p') : 0 : 0.0 : 1 : False : False : Binary
('d', 'q') : 0 : 0.0 : 1 : False : False : Binary
('d', 'r') : 0 : 0.0 : 1 : False : False : Binary
('d', 's') : 0 : 1.0 : 1 : False : False : Binary
a,b,c,d がそれぞれ p,q,r,sに対応しているのがわかる。
上のコードで以下のようにグラフの構造を変えてみる。
X = sort_order([("a", "b"), ("a", "c"), ("c", "d")])
Y = sort_order([("p", "q"), ("p", "r"), ("p", "s")])
すると解が得られず、この2つのグラフは同型でないことがわかる。
所感
非常に力ずくで、制約がかなりの数になるのが気になる。この最小の問題でも、成約数は50を超えている。ノード数が大きくなると破綻しそうだが、これでい良いんだろうか。