動画：グラフの分割問題（Python）
動画「QUBO++ ケーススタディ グラフの分割問題（Python）」の手引きです。 章ごとに、要点と動画に出てくるプログラム・出力、詳しい説明のページへのリンクをまとめています。 「▶ 動画の … から再生」をクリックすると、その場面から再生します。
この動画で扱うこと：
- 4 次元超立方体（頂点 16 個・辺 32 本）と、グラフを切ったときに切れる辺の数を考える 3 つの問題
- 3 つの問題に共通のカットの式
x[i] * ~x[j] + ~x[i] * x[j] - 最大カット（符号を反転して最小にする）と二分割幅（半分ずつという制約を
qbpp.cons()で加える） - カット幅：頂点の番号をドメインウォールで表し、変数の表の各列を切れ目での 2 つの分け方として、天井の整数変数
wで最大の切れ目を最小にする - 4 × 4 の盤のルークのグラフ（K4 × K4）で、辺の作り方だけを変えた同じプログラム
- 1 桁だけ違う番号の頂点を結ぶグラフでは番号順が最適であることと、K4 × K4 × K4 × K4（頂点 256 個）を
ABS3Solverで解く例
文章での説明はカット幅問題・最大カット問題・最小グラフ二分割問題のページにあります。動画（ケーススタディ）の一覧に戻る
4次元超立方体
- 頂点は 16 個で、0 から 15 の番号を 4 桁の 2 進数で書きます。ちょうど 1 ビットだけ違う頂点どうしを辺で結びます（例えば 0000 は 0001、0010、0100、1000 と結ばれます）。
- どの頂点も 4 つの頂点と結ばれているので、辺は 16 × 4 ÷ 2 = 32 本です。
- 最上位のビットが 0 の 8 個、1 の 8 個がそれぞれ 3 次元の立方体をつくり、2 つの立方体の対応する頂点どうしが 8 本の辺で結ばれています。
最大カット
- 頂点を 2 つのグループに分け、両端が違うグループにある辺を切れる辺、その数をカットと呼びます。カットが最大になる分け方を求めるのが最大カット問題です。
- 頂点 $i$ が一方のグループなら 1、もう一方なら 0 となる変数
x[i]を使います。辺 $(i, j)$ についてx[i] * ~x[j] + ~x[i] * x[j]は、両端が違うときだけ 1 です。すべての辺について足したものがカットの式cutです。 - ソルバーは最小にするので、
-cutを最小にします。 - 4 次元超立方体では 32 本すべてを切れます。辺の両端は 1 の個数が 1 つ違うので、1 の個数の偶数と奇数で分ければ、すべての辺が切れます。
プログラム（最大カットと二分割幅）
- 辺のリスト
edgesは、iのbビット目だけを反転した番号j = i ^ (1 << b)を使って作ります。i < jのときだけ加えるので、どの辺も 1 回ずつ入ります。 qbpp.varで 16 個のバイナリ変数の配列xを作り、カットの式cutを作ります。否定は変数の前の~です。- 最大カットは
f = -cutをsimplify_as_binary()で整理し、ExhaustiveSolverで解きます。変数が 16 個なので、すべての割り当てを調べられます。 - 二分割幅のプログラム（後半）は次の章で説明します。
import pyqbpp as qbpp
N = 16
edges = []
for i in range(N):
for b in range(4):
j = i ^ (1 << b) # i と 1 ビットだけ違う番号
if i < j:
edges.append((i, j))
M = len(edges)
x = qbpp.var("x", shape=N)
cut = 0
for i, j in edges:
cut += x[i] * ~x[j] + ~x[i] * x[j]
# 最大カット
f = -cut
f.simplify_as_binary()
sol = qbpp.ExhaustiveSolver(f).search()
print("max cut =", sol(cut))
print("S:", *[i for i in range(N) if sol(x[i]) == 1])
# 最小二分割
g = cut + (M + 1) * qbpp.cons(qbpp.sum(x) == N // 2)
g.simplify_as_binary()
sol = qbpp.ExhaustiveSolver(g).search()
print("bisection width =", sol(cut))
print("S:", *[i for i in range(N) if sol(x[i]) == 1])
出力：
max cut = 32
S: 1 2 4 7 8 11 13 14
bisection width = 8
S: 0 1 2 3 4 5 6 7
最大カットで表示されたグループは、2 進数の 1 の個数が奇数の頂点です。
詳しくは最大カット問題をご覧ください。
二分割幅
- 頂点を半分ずつの 2 つのグループに分けたときのカットの最小値を二分割幅と呼びます。
- 同じ
cutの式を使い、「xの和が 8」という条件をqbpp.cons()で制約にして、重みをかけてcutに足します（上のプログラムの後半）。 - 重みは辺の数に 1 を足した 33 です。カットは 32 以下なので、制約を破ると必ず損をします。
- 答えは 8 で、選ばれたのは 0〜7 の頂点です。最上位のビットで 2 つの立方体に分け、切れるのは立方体を結ぶ 8 本の辺だけです。
詳しくは最小グラフ二分割問題をご覧ください。
カット幅
- 頂点を 1 列に並べ、辺を頂点の上に配線します。同じレイヤーの辺どうしは重なってはいけません（端の頂点を共有するだけならかまいません）。
- 隣り合う頂点の間を切れ目と呼びます。1 つの切れ目を通る辺はすべて別のレイヤーにあるので、必要なレイヤーの数は、切れ目を通る辺の数の最大値です。これを並べ方のカット幅、すべての並べ方での最小値をグラフのカット幅と呼びます。
- 番号順に並べると、15 か所の切れ目を通る辺の数は次のとおりで、最大は 10 です。
4 6 8 8 10 10 10 8 10 10 10 8 8 6 4
- 真ん中の切れ目は頂点を 8 個ずつに分けます。二分割幅はあらゆる並べ方での真ん中の切れ目の最小値、カット幅は最大の切れ目の最小値なので、カット幅 ≥ 二分割幅です。
番号をドメインウォールで
- 頂点を 1 列に並べることは、左からの番号 0〜15（すべて異なる）を各頂点に割り当てることです。
- 頂点 $i$ の番号を 15 個のバイナリ変数のドメインウォール（1 が続いた後に 0 が続く形）で表し、壁の前の 1 の個数を番号とします。
x[i][k]は、頂点 $i$ の番号が $k$ より大きいときだけ 1 です。 - 変数
xは 16 行 15 列の 2 次元配列で、1 つの行が 1 つの頂点です。 - 1 つ目の制約：各行がドメインウォール（0 の次に 1 が来ない）。
- 2 つ目の制約：番号がすべて異なる。そのためには、列 $k$ の 1 の個数を $15-k$ にします。逆に、すべての列でこの個数なら、どの番号もちょうど 1 回ずつ使われます。
詳しくはカット幅問題の「ドメインウォールで番号を表す」をご覧ください。
切れ目のカットと天井
- 列 $k$ が 0 の頂点は番号が $k$ 以下、1 の頂点は $k$ より大きい頂点です。つまり列 $k$ は、$k$ 番目の切れ目で頂点を左右に分けています。
- そのため、切れ目を通る辺の数は最大カットと同じ式で数えられます。列 $k$ のカットの式を
cut[k]と書きます。 - 整数変数
wを用意し、すべての $k$ でcut[k]≤wという制約を課してwを最小にします。wはすべての切れ目のカットを上から押さえる天井で、下げていくと一番高い切れ目のところで止まります。 - 重みは、すべての制約で辺の数に 1 を足した 33 です。制約を破って得をできるのは
wが下がる分だけで、32 以下だからです。 - 変数は 240 個のバイナリ変数と 1 個の整数変数だけで、補助変数は使いません。
詳しくはカット幅問題の「列ごとのカット」「最大のカットを最小にする」をご覧ください。
プログラム（カット幅）
qbpp.varのshapeで 16 行 15 列の変数xを、integerで 0〜32 の整数変数wを作ります。- 1 つ目の制約は、スライスで列 1〜14 と列 0〜13 を取り出して引き、差が 0 以下とします。16 行 14 列の 224 本の制約が 1 行で書けます。
- 2 つ目の制約は、
qbpp.vector_sum()のaxis=0で並べた 15 個の列の和を、15 から 1 まで 1 ずつ減るリストと要素ごとに比べます。 - カットの式は最大カットと同じループですが、
x[i]は長さ 15 の行ベクトルなので、cutは 15 個の切れ目のカットが並んだ配列になります。cut - w <= 0の制約も 15 本まとめて書けます。 - 変数が 240 個あるので、
EasySolverに制限時間を与えて探索します。解から各行の 1 の個数を数えて頂点の番号を読み出し、番号の順に並べた頂点と各切れ目のカットを表示します。
import pyqbpp as qbpp
N = 16
edges = []
for i in range(N):
for b in range(4):
j = i ^ (1 << b) # i と 1 ビットだけ違う番号
if i < j:
edges.append((i, j))
M = len(edges)
x = qbpp.var("x", shape=(N, N - 1))
w = qbpp.var("w", integer=(0, M))
# 各行はドメインウォール
constraint = qbpp.cons(x[:, 1:] - x[:, :-1] <= 0)
# 番号がすべて異なる: 列 k の 1 は N-1-k 個
constraint += qbpp.cons(qbpp.vector_sum(x, axis=0) == list(range(N - 1, 0, -1)))
# 切れ目ごとのカット: x[i] は長さ N-1 の行ベクトル
cut = 0
for i, j in edges:
cut += x[i] * ~x[j] + ~x[i] * x[j]
constraint += qbpp.cons(cut - w <= 0)
f = w + (M + 1) * constraint
f.simplify_as_binary()
solver = qbpp.EasySolver(f)
sol = solver.search(time_limit=5)
print("cutwidth =", sol(w))
print("violated constraints =", f.cons(sol))
order = [0] * N
for i in range(N):
t = 0
for k in range(N - 1):
t += sol(x[i][k])
order[t] = i
print("order:", *order)
print("cuts:", *sol(cut))
実行例
出力の例：
cutwidth = 10
violated constraints = 0
order: 14 15 12 6 7 13 5 4 8 0 2 1 9 10 11 3
cuts: 4 6 8 10 10 10 10 8 10 10 10 10 8 6 4
- カット幅は 10、破った制約の数は 0 です。表示された順に頂点を並べると、10 レイヤーで配線できます。
- 並べ方は実行するたびに変わりますが、カット幅は 10 になります。
- 真ん中の切れ目を通るのは 8 本で、最大の 10 より少なくなっています。
詳しくはカット幅問題の「超立方体で確かめる」をご覧ください。
別のグラフで
- 4 × 4 の盤のマスを頂点とし、同じ行か同じ列のマスどうしを結びます（チェスのルークが 1 手で動けるマスどうし）。頂点は 16 個、辺は 48 本です。
- プログラムは、上の 2 つのプログラムの辺を作る部分だけを次のように変えます。マスの番号を 4 で割った商が行、余りが列です。
N = 16
edges = []
for i in range(N):
for j in range(i + 1, N):
if i // 4 == j // 4 or i % 4 == j % 4: # 同じ行か同じ列
edges.append((i, j))
M = len(edges)
最大カットと二分割幅の出力：
max cut = 32
S: 2 3 6 7 8 9 12 13
bisection width = 16
S: 0 1 2 3 4 5 6 7
カット幅の出力の例：
cutwidth = 18
violated constraints = 0
order: 8 12 15 13 14 9 10 11 5 4 7 3 2 0 6 1
cuts: 6 10 14 16 16 18 18 16 18 18 16 16 16 12 6
- 最大カットは 32 です。1 つの行の 4 つのマスは 6 本の辺で互いに結ばれ、2 個ずつに分けたときの 4 本が 1 つの行で切れる辺の最大です。行と列は合わせて 8 つなので、8 × 4 = 32 が上限で、出力された分け方はこれを達成しています。
- 二分割幅は 16 で、上の 2 行と下の 2 行に分けたときです。
- カット幅は 18 です。ここでも真ん中の切れ目は 16 本で、最大の 18 より少なくなります。
真ん中の切れ目が谷になる理由
- 番号順では、最上位の桁によって頂点がいくつかのまとまりに分かれ、真ん中の切れ目はちょうどまとまりの境目にあります。
- 境目を通るのはまとまりの間を結ぶ辺だけですが、境目から少しずれると、まとまりの内部も途中で切ることになるので、通る辺が多くなります。
- 二分割幅は真ん中だけを、カット幅はすべての切れ目を見るので、2 つは一致するとは限りません。
大きなグラフで
- 4 次元超立方体の頂点の番号は 2 進数の 4 桁、4 × 4 の盤のマスの番号は 4 進数の 2 桁（上の桁が行、下の桁が列）です。どちらも、1 桁だけ違う番号の頂点どうしを辺で結んだグラフです。
- この形のグラフでは、番号順の並べ方が最適であることが証明されています（K. Nakano, Linear Layouts of Generalized Hypercubes, WG 1993）。ソルバーが求めた 10 と 18 は番号順の並べ方のカット幅と同じなので、最適な値です。
- もっと大きな、4 進数の 4 桁のグラフ（K4 × K4 × K4 × K4）でも試します。頂点は 256 個、辺は 1536 本、変数の表は 256 行 255 列で、バイナリ変数は 65280 個です。
- プログラムは、辺を作る部分（4 進数の各桁を比べ、違う桁が 1 つだけなら辺を加える）と、ソルバーだけを変えます。ソルバーは GPU も使う
ABS3Solverです。
import pyqbpp as qbpp
def digits(v): # 4 進数の 4 桁
return [v // 4**t % 4 for t in range(4)]
N = 256
edges = []
for i in range(N):
for j in range(i + 1, N):
diff = sum(a != b for a, b in zip(digits(i), digits(j)))
if diff == 1: # 1 桁だけ違う
edges.append((i, j))
M = len(edges)
x = qbpp.var("x", shape=(N, N - 1))
w = qbpp.var("w", integer=(0, M))
# 各行はドメインウォール
constraint = qbpp.cons(x[:, 1:] - x[:, :-1] <= 0)
# 番号がすべて異なる: 列 k の 1 は N-1-k 個
constraint += qbpp.cons(qbpp.vector_sum(x, axis=0) == list(range(N - 1, 0, -1)))
# 切れ目ごとのカット: x[i] は長さ N-1 の行ベクトル
cut = 0
for i, j in edges:
cut += x[i] * ~x[j] + ~x[i] * x[j]
constraint += qbpp.cons(cut - w <= 0)
f = w + (M + 1) * constraint
f.simplify_as_binary()
solver = qbpp.ABS3Solver(f)
sol = solver.search(time_limit=10)
print("vertices =", N, " edges =", M)
print("cutwidth =", sol(w))
print("violated constraints =", f.cons(sol))
出力：
vertices = 256 edges = 1536
cutwidth = 306
violated constraints = 0
- 番号順の並べ方でも、切れ目を通る辺の数の最大は 306 なので、306 がこのグラフのカット幅です。
- 番号順の真ん中の切れ目は 256 本で、ここでも谷になっています。
まとめ
- 最大カット・二分割幅・カット幅の 3 つを、同じカットの式から作りました。
- 最大カットはカットの符号を反転して最小にし、二分割幅は半分ずつという制約を
qbpp.cons()で加えました。 - カット幅は頂点の番号をドメインウォールで表しました。変数の表の各列が切れ目での 2 つの分け方になり、天井の整数変数
wで最大の切れ目を最小にしました。 - 1 桁だけ違う番号の頂点を結ぶグラフでは番号順の並べ方が最適で、頂点 256 個の K4 × K4 × K4 × K4 でも、
ABS3Solverが最適なカット幅 306 を求めました。