最大次数制約を持つ最小スパニングツリー問題#
ここではJijZeptSolverとJijModelingを用いて、最大次数制約を持つ最小スパニングツリー問題を解く方法を説明します。 この問題は、Lucas, 2014, "Ising formulations of many NP problems" 8.1. Minimal Spanning Tree with a Maximal Degree Constraint で触れられています。
最大次数制約を持つ最小スパニングツリー問題とは#
最小スパニングツリー問題は、次のように定義されます。 各辺\((uv) \in E\)がそれぞれコスト\(c_{uv}\)を持つような、無向グラフ\(G=(V, E)\)が与えられたとします。 全ての頂点を含むようなツリーである、スパニングツリー \(T \subseteq G\)を構築することを考えます。 もし\(T\)が存在するとき、そのコスト
を最小化することを考えるのが、この最小スパニングツリー問題です。 そしてここに、\(T\)に含まれる頂点の次数が\(\Delta\)以下でなければならないという、最大次数制約を加えます。
例題#
簡単のため、ここでは次のような単純なグラフを見てみましょう。
左グラフにおいて、最大次数が3で制限されるような、最小スパニングツリーを求めてみましょう。
この問題に対する答えは、右のグラフのようになります。
このツリーの重みの総和は8で最小です。
そして最大次数制約も満たされていることがわかります。
数理モデル#
辺\(e=(uv)\)が最小スパニングツリー\(T\)に含まれるとき1、そうでないとき0となるようなバイナリ変数\(x_{uv}\)を導入します。
制約 1: 各頂点は最大で\(\Delta\)まで辺を持つことができる
グラフにおけるどの頂点の次数も、\(\Delta\)を超えることはできません。 ツリーでは少なくとも1つの辺を持たなければならないため、各頂点の次数に下界を考える必要があります。
\(x_{uv}, x_{vu}\)についての総和であることに注意が必要です。
制約 2: 選択される辺の数は\(|V|-1\)である
選択された辺の集合は、\(|V|\)個全ての頂点を含む木を形成しなければなりません。そのような全域木は、\(|V|-1\)本の辺を持ちます。
制約 3: 閉路の非存在
\(T\)がもし閉路\((v_1, v_2, \dots, v_k, v_1)\)を持つならば、その閉路は\(G\)の連結部分グラフを成しています。 したがって、\({v_1, v_2, \dots, v_k}\)を頂点とする部分グラフ\(S\)を考えると、\(S\)内で選ばれている辺の合計が\(k-1\)より大きくなります。
このような状況を排除するため、\(G\)の任意の連結部分グラフ\(S\)に対して、\(S\)内で選ばれている辺の合計が\(|S|-1\)以下でなければならないという制約を加えます。
ここで、\(E(S)\)は\(S\)内に含まれる辺の集合です。この制約と制約2を併せれば、選択された辺の集合が閉路を持たないことが保証されます。
一般に、\(|V|-1\)本の辺を持つ全域（＝頂点集合が\(V\)である）部分グラフ\(G'\)について、\(G'\)が閉路を持たないこと、\(G'\)が連結であることと、そして\(G'\)が全域木であることは全て同値になります。よって、制約2と制約3により、選択された辺の集合が全域木を形成することを保証できます。
目的関数
\(T\)の重みの総和を最小化する目的関数を設定します。
JijModelingによる定式化#
次に、JijModelingを用いて数理モデルを実装しましょう。 最初に、上述の数理モデルに必要な変数を定義します。
import jijmodeling as jm
# define problem
problem = jm.Problem("minimum spanning tree with a maximum degree constraint")
# define variables
V = problem.Natural('V') # number of vertices
E = problem.Graph('E', dtype=V) # set of edges
num_E = problem.NamedExpr('num_E', E.len_at(0))
C = problem.Float('C', dict_keys=E) # cost of each edge
D = problem.Natural('D') # maximum allowed degrees
# Each entry S[i] represents a connected subgraph and shall contain
# a list of (indices) of edges included in the subgraph, that is,
# if an edge E[k] is included in the subgraph represented by S[i], then k should be in S[i].
S = problem.Placeholder('S', dtype=num_E, shape=(None, num_E), jagged=True)
num_S = problem.NamedExpr('num_S', S.len_at(0))
# S_nodes[i] is the number of vertices in the subgraph represented by S[i].
S_nodes = problem.Natural('S_nodes', less_than=V, shape=(num_S,))
x = problem.BinaryVar('x', dict_keys=E)
制約と目的関数を実装しましょう。
problem += problem.Constraint(
'const1-1',
lambda v: jm.filter(lambda e: (e[0]==v)|(e[1]==v), E).map(lambda e: x[e]).sum() >= 1,
domain=V
)
problem += problem.Constraint(
'const1-2',
lambda v: jm.filter(lambda e: (e[0]==v)|(e[1]==v), E).map(lambda e: x[e]).sum() <= D,
domain=V
)
problem += problem.Constraint('const2', jm.sum(E, lambda e: x[e]) == V-1)
problem += problem.Constraint(
'const3',
lambda s: jm.sum(S[s], lambda e_i: x[E[e_i]]) <= S_nodes[s] - 1,
domain=num_S
)
problem += jm.sum(E, lambda e: C[e] * x[e])
Jupyter Notebook上で、実装した数理モデルを表示してみましょう。
problem
インスタンスの準備#
NetworkXを用いて、グラフを準備しましょう。 例として、ランダムな連結グラフを作成しましょう。
import itertools
import networkx as nx
import numpy as np
np.random.seed(seed=0)
# set the number of vertices
inst_V = 5
# set the number of degree
inst_D = 2
# set the probability of rewiring edges
p_rewire = 0.2
# set the number of nearest neighbors
k_neighbors = 4
# create a connected graph
inst_G = nx.connected_watts_strogatz_graph(inst_V, k_neighbors, p_rewire)
# add random costs to the edges
for u, v in inst_G.edges():
inst_G[u][v]['weight'] = np.random.randint(1, 10)
# create a list of edges
inst_E = [tuple(edge) for edge in inst_G.edges()]
# create a dictionary for the costs
inst_C = {(u, v): inst_G[u][v]['weight'] for u, v in inst_E}
# get subgraphs and their edges
subs_nodes = []
subs_edges = []
# get connected subgraphs that have at least 2 nodes and have less nodes than the original graph
for num_nodes in range(2, inst_G.number_of_nodes()):
for comb in (inst_G.subgraph(selected_nodes) for selected_nodes in itertools.combinations(inst_G, num_nodes)):
if nx.is_connected(comb):
edges_indices = [inst_E.index((u, v)) for u, v in comb.edges()]
subs_edges.append(edges_indices)
subs_nodes.append(comb.nodes)
inst_S = subs_edges
inst_S_nodes = [len(subgraph_nodes) for subgraph_nodes in subs_nodes]
instance_data = {'V': inst_V, 'D': inst_D, 'E': inst_E, 'C': inst_C, 'S': inst_S, 'S_nodes': inst_S_nodes}
ここで作られたグラフを表示してみます。
import matplotlib.pyplot as plt
pos = nx.spring_layout(inst_G)
nx.draw_networkx(inst_G, pos, with_labels=True)
edge_labels = nx.get_edge_attributes(inst_G, 'weight')
nx.draw_networkx_edge_labels(inst_G, pos, edge_labels=edge_labels, font_color='red')
plt.show()
Solve by JijZeptSolver#
この問題をjijzept_solverを用いて解きましょう。
import jijzept_solver
instance = problem.eval(instance_data)
solution = jijzept_solver.solve(instance, solve_limit_sec=1.0)
解の可視化#
得られた解を可視化してみましょう。
# get the indices of x == 1
df = solution.decision_variables_df
x_indices = df[df["value"]==1]["subscripts"].to_list()
# draw figure
plt.subplot(121)
plt.axis('off')
nx.draw(inst_G, pos, with_labels=True, font_weight='bold')
nx.draw_networkx_edge_labels(inst_G, pos, edge_labels=edge_labels, font_color='red')
plt.subplot(122)
plt.axis('off')
inst_H = inst_G.copy()
for e in df[(df["name"]=="x") & (df["value"] == 0)]["subscripts"]:
inst_H.remove_edge(*e)
nx.draw(inst_H, pos, with_labels=True, font_weight='bold')
nx.draw_networkx_edge_labels(inst_H, pos, edge_labels=nx.get_edge_attributes(inst_H, 'weight'), font_color='red')
{(0, 4): Text(0.004709307600212931, 0.0906234383694593, '1'),
(0, 2): Text(0.49403114046146035, 0.5678539174374728, '4'),
(1, 3): Text(-0.10110844399733399, -0.15848411239151816, '4'),
(2, 3): Text(0.3237174278353945, -0.36524810117169304, '3')}
期待通り、最大次数制約を満たす最小スパニングツリーを得ることができました。
ここで紹介した問題の高度な発展として、シュタイナーツリー問題というものがあります。 この問題では、辺のコスト\(c_{uv}\)と頂点の部分集合\(U \subset V\)が与えられた時、 \(U\)に含まれる全ての頂点を含むような最小コストの木を求めます。
シュタイナーツリー問題を解くには最小スパニングツリー問題と同じハミルトニアンを用いれば良いですが、部分集合\(U\)に含まれない頂点\(v\)を木に含めるかどうかを表す追加のバイナリ変数\(y_v\)を用いる必要があります。