グラフ2分割問題
グラフ2分割問題(graph bipartition problem)は、以下のように定義される。
点数 $n=|V|$ が偶数である無向グラフ $G=(V,E)$ が与えられたとき， 点集合 $V$ の等分割(uniform partition, eqipartition) $(L,R)$ とは， $L \cap R = \emptyset$, $L \cup R =V$, $|L|=|R|=n/2$ を満たす点の部分集合の対である．
$L$ と $R$ をまたいでいる枝 （正確には $i \in L, j \in R$ もしくは $i \in R, j \in L$ を満たす枝 $(i,j)$） の本数を最小にする等分割 $(L,R)$ を求める問題がグラフ2分割問題である．
一般には $k (>2)$ 分割も考えられるが，ここでは $k=2$ に絞って議論する．
問題を明確化するために，グラフ分割問題を整数最適化問題として定式化しておく． 無向グラフ $G=(V,E)$ に対して，$L \cap R =\emptyset$（共通部分がない）， $L \cup R =V$（合わせると点集合全体になる）を満たす非順序対 $(L,R)$ を2分割(bipartition)と呼ぶ． 分割 $(L,R)$ において，$L$ は左側，$R$ は右側を表すが，これらは逆にしても同じ分割であるので，非順序対と呼ばれる．
点 $i$ が，分割 $(L,R)$ の $L$ 側に含まれているとき $1$， それ以外の（$R$ 側に含まれている）とき $0$ の $0$-$1$ 変数 $x_i$ を導入する． このとき，等分割であるためには，$x_i$ の合計が $n/2$ である必要がある． 枝 $(i,j)$ が $L$ と $R$ をまたいでいるときには，$x_i (1-x_j)$ もしくは $(1-x_i)x_j$ が $1$ になることから，以下の定式化を得る．
$$ \begin{array}{l l l} minimize & \sum_{(i,j) \in E} x_i (1-x_j)+(1-x_i)x_j & \\ s.t. & \sum_{i \in V} x_i = n/2 & \\ & x_i \in \{0,1\} & \forall i \in V \end{array} $$数理最適化ソルバーの多くは，上のように凸でない2次の項を目的関数に含んだ最小化問題には対応していない． したがって，一般的な（混合）整数最適化ソルバーで求解するためには，2次の項を線形関数に変形してから解く必要がある．
枝 $(i,j)$ が $L$ と $R$ をまたいでいるとき $1$，それ以外のとき $0$ になる $0$-$1$ 変数 $y_{ij}$ を導入する． すると，上の2次整数最適化問題は，以下の線形整数最適化に帰着される．
$$ \begin{array}{l l l} minimize & \sum_{(i,j) \in E} y_{ij} & \\ s.t. & \sum_{i \in V} x_i = n/2 & \\ & x_i -x_j \leq y_{ij} & \forall (i,j) \in E \\ & x_j -x_i \leq y_{ij} & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \\ & y_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array} $$最初の制約は等分割を規定する． 2番目の制約は，$i \in L$ で $j \not\in L$ のとき $y_{ij}=1$ になることを規定する． 3番目の制約は，$j \in L$ で $i \not\in L$ のとき $y_{ij}=1$ になることを規定する．
上の定式化を行う関数を以下に示す．
V, E = make_data(10, 0.4)
model = gpp(V, E)
model.optimize()
print("Opt.value=", model.ObjVal)
x, y = model.__data
L = [i for i in V if x[i].X >= 0.5]
R = [i for i in V if x[i].X < 0.5]
G = nx.Graph()
G.add_edges_from(E)
pos = nx.layout.spring_layout(G)
nx.draw(G, pos=pos, with_labels=True, node_size=500, node_color="yellow", nodelist=L)
edgelist = [(i, j) for (i, j) in E if (i in L and j in R) or (i in R and j in L)]
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=500,
nodelist=R,
node_color="Orange",
edgelist=edgelist,
edge_color="red",
)
# ampl.option['solver'] = 'highs'
# ampl.read("../ampl/gpp.mod")
# ampl.set['V'] = list(V)
# ampl.set['E'] = list(E)
# ampl.solve()
# z = ampl.obj['z']
# print("Optimum:", z.value())
# x = ampl.var['x']
# L = [i for i in V if x[i].value() >.5]
# R = [i for i in V if x[i].value() <.5]
# print("L =", L)
# print("R =", R)
# y = ampl.var['y']
# across = [(i,j) for (i,j) in E if y[i,j].value() > .5]
# print("Across edges:", across)
rndseed = 2
random.seed(rndseed)
n = 100
pos = {i: (random.random(), random.random()) for i in range(n)}
G = nx.random_geometric_graph(n, 0.2, pos=pos)
nodes = G.nodes()
adj = [set([]) for i in nodes]
for (i, j) in G.edges():
adj[i].add(j)
adj[j].add(i)
max_iterations = 1000
tabulen = len(nodes) / 2
sol = construct(nodes)
z, _, s, d = evaluate(nodes, adj, sol)
print("initial partition: z =", z)
print(sol)
print()
print("starting tabu search,", max_iterations, "iterations, tabu length =", tabulen)
sol, cost = tabu_search(nodes, adj, sol, max_iterations, tabulen)
print("final solution: z=", cost)
print(sol)
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
(i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
nodelist=R,
node_color="Orange",
edgelist=edgelist,
edge_color="red",
)
アニーリング法
ランダム性を用いて局所的最適解からの脱出を試みるメタヒューリスティクスの代表例として アニーリング法(simulated annealing method）がある． この解法の土台は，1953年に Metropolis--Rosenbluth--Rosenbluth--Teller--Tellerによって築かれたものであり，何度も再発見されている． そのため，様ざなま別名をもつ．一部を紹介すると， モンテカルロ焼なまし法（Monte Carlo annealing）, 確率的丘登り法（probabilistic hill climbing）, 統計的冷却法（statistical cooling）, 確率的緩和法（stochastic relaxation）, Metropolisアルゴリズム（Metropolis algorithm） などがある．
アニーリング法の特徴は，温度とよばれるパラメータを徐々に小さくしていくことによって，改悪の確率を $0$ に近づけていくことである．． 以下では，分割の非均衡度 $|L|-|R|$ をペナルティとしたアニーリング法を示す．
initprob = 0.4 # initial acceptance probability
sizefactor = 16 # for det. # tentatives at each temp.
tempfactor = 0.95 # cooling ratio
freezelim = 5 # max number of iterations with less that minpercent acceptances
minpercent = 0.02 # fraction of accepted moves for being not frozen
alpha = 3.0 # imballance factor
N = len(nodes) # neighborhood size
L = N * sizefactor # number of tentatives at current temperature
print("starting simulated annealing, parameters:")
print(" initial acceptance probability", initprob)
print(" cooling ratio", tempfactor)
print(" # tentatives at each temp.", L)
print(" percentage of accepted moves for being not frozen", minpercent)
print(" max # of it.with less that minpercent acceptances", freezelim)
print(" imballance factor", alpha)
print()
sol = construct(nodes) # starting solution
z, bal, s, d = evaluate(nodes, adj, sol, alpha)
print("initial partition: z =", z)
print(sol)
print()
LOG = False
sol, z = annealing(
nodes, adj, sol, initprob, L, tempfactor, freezelim, minpercent, alpha, report=False
)
zp, bal, sp, dp = evaluate(nodes, adj, sol, alpha)
assert abs(zp - z) < 1.0e-9 # floating point approx.equality
print()
print("final solution: z=", z)
print(sol)
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
(i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
nodelist=R,
node_color="Orange",
edgelist=edgelist,
edge_color="red",
)
sol = construct(nodes)
print("initial partition: z =", z)
print(sol)
print()
print("starting tabu search,", max_iterations, "iterations, tabu length =", tabulen)
sol, cost = ts_intens_divers(nodes, adj, sol, max_iterations, tabulen, report=False)
print("final solution: z=", cost)
print(sol)
# G.nodes[v]["color"] = sol[v]
# fig = gts.to_plotly_fig(G,node_size=30,
# line_width=0.1,
# line_color='blue',
# text_size=10, pos=pos)
# plotly.offline.plot(fig);
L = [i for i in nodes if sol[i] == 1]
R = [i for i in nodes if sol[i] == 0]
edgelist = [
(i, j) for (i, j) in G.edges() if (i in L and j in R) or (i in R and j in L)
]
nx.draw(G, pos=pos, with_labels=False, node_size=100, node_color="yellow", nodelist=L)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
nodelist=R,
node_color="Orange",
edgelist=edgelist,
edge_color="red",
)
rndseed = 2
random.seed(rndseed)
n = 10000
pos = {i: (random.random(), random.random()) for i in range(n)}
G = nx.random_geometric_graph(n, 0.05, pos=pos)
#枝に重みを付与
for u, v in G.edges():
G[u][v]["weight"] = random.randint(1,1)
#点に重みを付与
for u in G.nodes():
G.nodes[u]["weight"] = random.randint(1,100)
nx.draw(G, pos=pos, node_size=0.3)
npartition = 5 #分割数
ufactor = 30 #不均衡パラメータ
sol = metis(G, npartition, ufactor)
num_nodes = np.zeros(npartition, int)
for i in sol:
num_nodes[sol[i]] += G.nodes[i]["weight"]
color_list = [sol[i] for i in G.nodes()]
nx.draw_networkx_nodes(G, pos, node_color=color_list, node_size=5)
# 枝の描画
edge_colors = []
edge_widths = []
obj = 0
for u, v in G.edges():
if sol[u] != sol[v]: # 両端ノードの色が異なる場合
obj += G[u][v]["weight"]
edge_colors.append("r") # 赤色で描画
edge_widths.append(1) # 太さを2に設定
else:
edge_colors.append("k") # 黒色で描画
edge_widths.append(0.1) # 太さを1に設定
nx.draw_networkx_edges(G, pos, edge_color=edge_colors, width=edge_widths)
print(obj, num_nodes)
線形定式化
グラフ分割問題と同様に， 点 $i$ が，分割 $(L,R)$ の $L$ 側に含まれているとき $1$， それ以外の（$R$ 側に含まれている）とき $0$ の $0$-$1$ 変数 $x_i$ を導入する． 枝 $(i,j)$ が $L$ と $R$ をまたいでいるとき $1$，それ以外のとき $0$ になる $0$-$1$ 変数 $y_{ij}$ を導入する．
2分割を表す制約のかわりに， 点 $i,j$ が同じ側に含まれているときに $y_{ij}$ を $0$ にするための制約を追加する．
$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} w_{ij} y_{ij} & \\ s.t. & x_i + x_j \geq y_{ij} & \forall (i,j) \in E \\ & 2- x_i - x_j \geq y_{ij} & \forall (i,j) \in E \\ & x_i -x_j \leq y_{ij} & \forall (i,j) \in E \\ & x_j -x_i \leq y_{ij} & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \\ & y_{ij} \in \{0,1\} & \forall (i,j) \in E \end{array} $$def maxcut(V, E):
"""maxcut -- model for the graph maxcut problem
Parameters:
- V: set/list of nodes in the graph
- E: set/list of edges in the graph
Returns a model, ready to be solved.
"""
model = Model("maxcut")
x = {}
y = {}
for i in V:
x[i] = model.addVar(vtype="B", name=f"x({i})")
for (i, j) in E:
y[i, j] = model.addVar(vtype="B", name=f"y({i},{j})")
model.update()
for (i, j) in E:
model.addConstr(x[i] + x[j] >= y[i, j], f"Edge({i},{j})")
model.addConstr(2 - x[j] - x[i] >= y[i, j], f"Edge({j},{i})")
model.addConstr(x[i] - x[j] <= y[i, j], f"EdgeLB1({i},{j})")
model.addConstr(x[j] - x[i] <= y[i, j], f"EdgeLB2({j},{i})")
model.setObjective(quicksum(y[i, j] for (i, j) in E), GRB.MAXIMIZE)
model.update()
model.__data = x
return model
2次錐最適化による定式化
ここでは，以下の論文による2次錐最適化による強い定式化を示す．
A NEW SECOND-ORDER CONE PROGRAMMING RELAXATION FOR MAX-CUT PROBLEMS, Masakazu Muramatsu Tsunehiro Suzuki, Journal of the Operations Research, 2003, Vol. 46, No. 2, 164-177
$x_j$ は上の定式化と同じであるが，以下の2つの $0$-$1$ 変数を導入する．
- $s_{ij}$: $i$ と $j$ が同じ分割に含まれるとき $1$
- $z_{ij}$: $i$ と $j$ が異なる分割に含まれるとき $1$
def maxcut_soco(V, E):
"""maxcut_soco -- model for the graph maxcut problem
Parameters:
- V: set/list of nodes in the graph
- E: set/list of edges in the graph
Returns a model, ready to be solved.
"""
model = Model("max cut -- scop")
x, s, z = {}, {}, {}
for i in V:
x[i] = model.addVar(vtype="B", name=f"x({i})")
for (i, j) in E:
s[i, j] = model.addVar(vtype="C", name=f"s({i},{j})")
z[i, j] = model.addVar(vtype="C", name=f"z({i},{j})")
#model.update()
for (i, j) in E:
model.addConstr((x[i] + x[j] - 1) * (x[i] + x[j] - 1) <= s[i, j], f"S({i},{j})")
model.addConstr((x[j] - x[i]) * (x[j] - x[i]) <= z[i, j], f"Z({i},{j})")
model.addConstr(s[i, j] + z[i, j] == 1, f"P({i},{j})")
model.setObjective(quicksum(z[i, j] for (i, j) in E), GRB.MAXIMIZE)
#model.update()
#model.__data = x
return model
def make_data(n, prob):
"""make_data: prepare data for a random graph
Parameters:
- n: number of vertices
- prob: probability of existence of an edge, for each pair of vertices
Returns a tuple with a list of vertices and a list edges.
"""
V = range(1, n + 1)
E = [(i, j) for i in V for j in V if i < j and random.random() < prob]
return V, E
#GRB = MDO
from gurobipy import Model, quicksum, GRB
random.seed(3)
V, E = make_data(100, 0.1)
model = maxcut(V, E)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
print("Opt.value=", model.ObjVal)
obj_lin = model.ObjVal
x = model.__data
print([i for i in V if x[i].X >= 0.5])
print([i for i in V if x[i].X < 0.5])
else:
pass
#GRB = MDO
from gurobipy import Model, quicksum, GRB
model = maxcut_soco(V, E)
model.optimize()
status = model.Status
if status == GRB.Status.OPTIMAL:
print("Opt.value=", model.ObjVal)
x = model.__data
print([i for i in V if x[i].X >= 0.5])
print([i for i in V if x[i].X < 0.5])
# assert model.ObjVal == obj_lin
else:
pass
2次無制約2値最適化による定式化
Gurobiの最近のバージョンでは，2次無制約2値最適化 (Quadratic unconstrained binary optimization: QUBO)に対する専用のヒューリスティクスを組み込んでいる．
最大カット問題は，QUBOとして以下のように定式化できる． QUBOでは，2値変数は $-1$ か $1$ の値ととるものと仮定する． 分割 $(L,R)$ の $L$ 側に含まれているとき $1$， それ以外の（$R$ 側に含まれている）とき $-1$ の2値変数 $x_i$ を用いる．
$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} w_{ij} (1-x_{i}x_{j}) & \\ & x_i \in \{-1,1\} & \forall i \in V \end{array} $$これを，以下の変換を用いて $0$-$1$ 変数 $\hat{x}_i$ に変換する．
$$ \hat{x}_i = \frac{x_i+1}{2} $$これをQUBOによる代入すると，2次無制約の$0$-$1$ 整数最適化は，以下のようになる．
$$ \begin{array}{l l l} maximize & \sum_{(i,j) \in E} \frac{1}{2} w_{ij} (-4 \hat{x}_{i} \hat{x}_{j}+ 2\hat{x}_{i} + 2 \hat{x}_{j}) & \\ & \hat{x}_i \in \{0,1\} & \forall i \in V \end{array} $$以下のサイトからベンチマーク問題例を入手できるので， ベンチマーク問題例で実験する．
from gurobipy import Model, quicksum, GRB
folder = "../data/maxcut/"
gurobi = []
for iter in range(1,2): #ここで問題番号を指定する(全部で54問)
f = open(folder + f"g{iter}.rud")
data = f.readlines()
f.close()
n, m = list(map(int, data[0].split()))
G = nx.Graph()
edges = []
for row in data[1:]:
i, j, w = list(map(int, row.split()))
edges.append((i, j, w))
G.add_weighted_edges_from(edges)
model = Model()
x = {}
for i in G.nodes():
x[i] = model.addVar(name=f"x[{i}]",vtype="B")
model.setObjective( quicksum(1/2* G[i][j]["weight"]*(-4*x[i]*x[j]+2*x[i]+2*x[j]) for (i, j) in G.edges()), GRB.MAXIMIZE)
model.Params.TimeLimit = 10
model.Params.OutputFlag=1
model.optimize()
print(iter, model.ObjVal)
gurobi.append(model.ObjVal)
制約最適化ソルバー SCOP による求解
上と同じ問題例に対して，制約最適化ソルバー SCOP （付録1参照）で求解を試みる．
変数は $X_i (i \in V)$ で値集合は $\{0,1\}$ とする． 値変数は $x_{i0}, x_{i1}$ となる．
ベンチマーク問題例には，枝の重み $w_{ij}, (i,j) \in E$が $1$ のものと $-1$ のものがある．
枝の重みが $1$ の枝に対しては， 重み（逸脱ペナルティ）が $1$ の考慮制約として，以下の2次制約を定義する．
$$ x_{i0} x_{j1} + x_{i1} x_{j0} \geq 1 \ \ \ \forall (i,j) \in E, w_{ij}=1 $$枝 $(i,j)$ の両端点が異なる集合に含まれるとき左辺は $1$ となり，制約は満たされる． それ以外のとき左辺は $0$ となり，制約は $1$ だけ破られる．
枝の重みが $-1$ の枝に対しては， 重み（逸脱ペナルティ）が $1$ の考慮制約として，以下の2次制約を定義する．
$$ x_{i0} x_{j1} + x_{i1} x_{j0} \leq 0 \ \ \ \forall (i,j) \in E, w_{ij}=-1 $$枝 $(i,j)$ の両端点が異なる集合に含まれるとき左辺は $1$ となり，制約は $1$ だけ破られる．
枝の本数 $|E|$ から逸脱の総数を減じたものが，目的関数になる．
SCOPモジュールからインポートする．
from scop import *
folder = "../data/maxcut/"
for iter in range(1, 2): #ここで問題番号を指定する(全部で54問)
f = open(folder + f"g{iter}.rud")
data = f.readlines()
f.close()
n, m = list(map(int, data[0].split()))
G = nx.Graph()
edges = []
for row in data[1:]:
i, j, w = list(map(int, row.split()))
edges.append((i, j, w))
G.add_weighted_edges_from(edges)
model = Model()
x = {}
for i in G.nodes():
x[i] = model.addVariable(name=f"x[{i}]", domain=[0, 1])
for (i, j) in G.edges():
if G[i][j]["weight"] == 1:
q = Quadratic(name=f"edge_{i}_{j}", weight=1, rhs=1, direction=">=")
q.addTerms(
coeffs=[1, 1],
vars=[x[i], x[i]],
values=[0, 1],
vars2=[x[j], x[j]],
values2=[1, 0],
)
model.addConstraint(q)
elif G[i][j]["weight"] == -1:
q = Quadratic(name=f"edge_{i}_{j}", weight=1, rhs=0, direction="<=")
q.addTerms(
coeffs=[1, 1],
vars=[x[i], x[i]],
values=[0, 1],
vars2=[x[j], x[j]],
values2=[1, 0],
)
model.addConstraint(q)
else:
print("no such edge! ")
break
model.Params.TimeLimit = 100
model.Params.RandomSeed = 123
sol, violated = model.optimize()
val = 0
for (i, j) in G.edges():
if x[i].value != x[j].value:
if G[i][j]["weight"] == 1:
val += 1
elif G[i][j]["weight"] == -1:
val -= 1
print(iter, n, m, val)
SCOPとGurobi （ver. 10) ともに100秒で打ち切ったときの実験結果を以下に示す．
folder = "../data/maxcut/"
df = pd.read_csv(folder + "maxcut_bestvalues.csv", index_col=0)
df["gap"] = df["Best Known"] / df["scop(100s)"] * 100
df
#df["gurobi"] = gurobi
#df.to_csv(folder + "maxcut_bestvalues.csv")