最大クリーク問題(最大安定集合問題)とその周辺

準備

ここでは,付録2で準備したグラフに関する基本操作を集めたモジュール graphtools.py を読み込んでいる. 環境によって,モジュールファイルの場所は適宜変更されたい.

最大クリーク問題と最大安定集合問題

最大安定集合問題(maximum stable set problem)は,以下のように定義される問題である.

点数 $n$ の無向グラフ $G=(V,E)$ が与えられたとき,点の部分集合 $S (\subseteq V)$ は, すべての $S$ 内の点の間に枝がないとき安定集合(stable set)とよばれる. 最大安定集合問題とは,集合に含まれる要素数(位数)$|S|$ が最大になる安定集合 $S$ を求める問題である.

この問題のグラフの補グラフ(枝の有無を反転させたグラフ)を考えると, 以下に定義される最大クリーク問題(maximum clique problem)になる.

無向グラフ $G=(V,E)$ が与えられたとき,点の部分集合 $C (\subseteq V)$は, $C$ によって導かれた誘導部分グラフが完全グラフ(complete graph)になるときクリーク(clique)とよばれる (完全グラフとは,すべての点の間に枝があるグラフである). 最大クリーク問題とは,位数 $|C|$ が最大になるクリーク $C$ を求める問題である.

これらの2つの問題は(お互いに簡単な変換によって帰着されるという意味で)同値である.

点 $i$ が安定集合 $S$ に含まれるとき $1$,それ以外のとき $0$ の $0$-$1$ 変数を用いると, 最大安定集合問題は,以下のように定式化できる.

$$ \begin{array}{l l l} maximize & \sum_{i \in V} x_i & \\ s.t. & x_i +x_j \leq 1 & \forall (i,j) \in E \\ & x_i \in \{0,1\} & \forall i \in V \end{array} $$

極大クリークの列挙

networkXのfind_cliquesを使うことによってグラフの極大クリークを列挙できる. パスの列挙で紹介した Graphillionは,位数を固定して列挙できるが,遅いので使うべきではない.

また,グラフが疎でないと極大クリーク数は指数関数的に増大するので,注意が必要である.

G = nx.fast_gnp_random_graph(100, 0.5)
count = 0
for i in nx.find_cliques(G):
    count += 1
count
16743
count = 0
for i in nx.find_cliques(G):
    print(i)
    count += 1
    if count >= 10:
        break
[1, 65, 32, 75, 73, 3]
[1, 65, 32, 75, 73, 36, 38]
[1, 65, 32, 75, 91, 3]
[1, 65, 32, 75, 20]
[1, 65, 32, 83, 80, 38]
[1, 65, 32, 98, 26, 91]
[1, 65, 32, 98, 26, 84, 73]
[1, 65, 32, 98, 26, 38, 80]
[1, 65, 32, 98, 26, 38, 73]
[1, 65, 32, 98, 36, 73, 84]

近似解法

networkXに $O(|V|/(log|V|)^2)$ の近似精度保証をもった近似解法があるが,大きな問題例だと遅い. また,実験結果も思わしくないので,使うべきではない. 往々にして,性能保証をもった近似解法の性能は悪い.

100点のランダムなグラフに対して近似解法を適用する.

n = 100
random.seed(1)
pos = {i: (random.random(), random.random()) for i in range(n)}
G = nx.random_geometric_graph(n, 0.3, pos=pos, seed=1)
from networkx.algorithms import approximation
S = approximation.max_clique(G)
print(len(S))
13
nx.draw(G, pos=pos, node_size=100)
nx.draw(
    G,
    pos=pos,
    nodelist=list(S),
    edgelist=[(i, j) for i in S for j in S if i < j],
    node_color="red",
    edge_color="blue",
)

タブーサーチ

最大安定集合を求めるためのメタヒューリスティクスとして,タブーサーチを設計する.

以下のメタヒューリスティクスの詳細については,拙著「メタヒューリスティクスの数理」(共立出版)を参照されたい.

evaluate[source]

evaluate(nodes, adj, sol)

Evaluate a solution.

Determines:

  • the cardinality of a solution, i.e., the number of nodes in the stable set;
  • the number of conflicts in the solution (pairs of nodes for which there is an edge);
  • b[i] - number of nodes adjacent to i in the stable set;

construct[source]

construct(nodes, adj)

A simple construction method.

The solution is represented by a set, which includes the vertices in the stable set.

This function constructs a maximal stable set.

find_add[source]

find_add(nodes, adj, sol, b, tabu, tabulen, iteration)

Find the best non-tabu vertex for adding into 'sol' (the stable set)

find_drop[source]

find_drop(nodes, adj, sol, b, tabu, tabulen, iteration)

Find the best non-tabu vertex for removing from 'sol' (the stable set)

move_in[source]

move_in(nodes, adj, sol, b, tabu, tabuIN, tabuOUT, iteration)

Determine and execute the best non-tabu insertion into the solution.

move_out[source]

move_out(nodes, adj, sol, b, tabu, tabuIN, tabuOUT, iteration)

Determine and execute the best non-tabu removal from the solution.

tabu_search(nodes, adj, sol, max_iter, tabulen, report=None)

Execute a tabu search run.

Gbar = nx.complement(G)
nodes, edges = Gbar.nodes(), Gbar.edges()
adj = gts.adjacent(nodes, edges)
sol = construct(nodes, adj)

max_iter = 10000
tabulen = len(nodes) / 10
bestsol, bestcard = tabu_search(nodes, adj, sol, max_iter, tabulen, report=None)

nx.draw(G, pos=pos, node_size=100)
nx.draw(
    G,
    pos=pos,
    nodelist=list(S),
    edgelist=[(i, j) for i in S for j in S if i < j],
    node_color="red",
    edge_color="blue",
)
print(len(bestsol))
15

集中化・多様化を入れたタブーサーチ

上のタブーサーチに集中化と多様化を加味して改善する.

diversify[source]

diversify(nodes, adj, v)

Find a maximal stable set, starting from node 'v'

ts_intens_divers[source]

ts_intens_divers(nodes, adj, sol, max_iter, tabulen, report)

Execute a tabu search run using intensification/diversification.

sol = construct(nodes, adj)
max_iter = 10000
tabulen = len(nodes) / 10
bestsol, bestcard = ts_intens_divers(nodes, adj, sol, max_iter, tabulen, report=False)

nx.draw(G, pos=pos, node_size=100)
nx.draw(
    G,
    pos=pos,
    nodelist=list(S),
    edgelist=[(i, j) for i in S for j in S if i < j],
    node_color="red",
    edge_color="blue",
)
print(len(bestsol))
15

平坦探索法

最大安定集合問題(最大クリーク問題)は,同じ目的関数をもつ解がたくさんあるという特徴をもつ.これを平坦部とよぶ. この平坦部からの脱出に工夫を入れたメタヒューリスティクスを設計する.

possible_add[source]

possible_add(rmn, b)

Check which of the nodes in set 'rmn' can be added to the current stable set (i.e., which have no edges to nodes in the stable set, and thus have b[i] = 0).

one_edge[source]

one_edge(rmn, b)

Check which of the nodes in 'rmn' cause one conflict if added to the current stable set (i.e., which have exactly one edge to nodes in the stable set, and thus have b[i] = 1).

add_node[source]

add_node(i, adj, sol, rmn, b)

Move node 'i' from 'rmn' into 'sol', and update 'b' accordingly.

expand_rand[source]

expand_rand(add, sol, rmn, b, adj, maxiter, dummy=None)

Expand the current stable set ('sol') from randomly selected valid nodes. Use nodes in 'add' for the expansion; 'add' must be the subset of unselected nodes that may still be added to the stable set.

expand_stat_deg[source]

expand_stat_deg(add, sol, rmn, b, adj, maxiter, degree)

Expand the current stable set ('sol'), selecting nodes with maximal degree on the initial graph. Use nodes in 'add' for the expansion; 'add' must be the subset of unselected nodes that may still be added to the stable set.

expand_dyn_deg[source]

expand_dyn_deg(add, sol, rmn, b, adj, maxiter, dummy=None)

Expand the current stable set ('sol'), selecting nodes with maximal degree on 'add' (i.e., the subset of nodes that can still be added to the stable set).

iterated_expansion[source]

iterated_expansion(nodes, adj, expand_fn, niterations, report=None)

Do repeated expansions until reaching 'niterations'. Expansion is done using the function 'expand_fn' coming as a parameter.

Parameters:

  • nodes, adj -- graph information
  • expand_fn -- function to be used for the expansion
  • niterations -- maximum number of iterations

node_replace[source]

node_replace(v, sol, rmn, b, adj)

Node 'v' has been inserted in 'sol' and created one conflict. Remove the conflicting node (the one in 'sol' adjacent to 'v'), update 'b', and the check through which nodes expansion is possible.

plateau[source]

plateau(sol, rmn, b, adj, maxiter)

Check nodes that create one conflict if inserted in the stable set; tentatively add them, and remove the conflict created. Exit whenever the stable set can be expanded (and return the subset of nodes usable for that).

multistart_local_search(nodes, adj, expand_fn, niterations, length, report=None)

Plateau search, using 'expand_fn' for the expansion.

Parameters:

  • nodes, adj -- graph information
  • expand_fn -- function to be used for the expansion
  • niterations -- maximum number of iterations
  • length -- number of searches to do attempt on the each plateau

expand_through[source]

expand_through(add, sol, rmn, expand_fn, b, adj, maxiter, degree)

Expand the current stable set ('sol'). Initially use nodes in 'add' for the expansion; then, try with all the unselected nodes (those in 'rmn').

ltm_search(nodes, adj, expand_fn, niterations, length, report=None)

Plateau search including long-term memory, using 'expand_fn' for the expansion.

Long-term memory 'ltm' keeps the number of times each node was used after successful expansions.

Parameters:

  • nodes, adj -- graph information
  • expand_fn -- function to be used for the expansion
  • niterations -- maximum number of iterations
  • length -- number of searches to do attempt on the each plateau

rm_node[source]

rm_node(i, adj, sol, rmn, b)

Move node 'i' from 'sol' into 'rmn', and update 'b' accordingly.

hybrid[source]

hybrid(nodes, adj, niterations, length, report=None)

Plateau search, using a hybrid approach.

The solution is partially kept after each plateau search: a node from the unselected 'rmn' set is chosen and added to the stable set, and then the conflicting nodes are removed.

Plateau expansions are done through a randomly selected method, from random expansion, dynamic-degree-based expansion, or static-degree expansion.

Long-term memory 'ltm' keeps the number of times each node was used after successful expansions.

sol = construct(nodes, adj)
max_iterations = 10000
tabulength = len(nodes) / 10
sol, rmn, card = hybrid(nodes, adj, max_iterations, tabulength, False)
nx.draw(G, pos=pos, node_size=100)
nx.draw(
    G,
    pos=pos,
    nodelist=list(S),
    edgelist=[(i, j) for i in S for j in S if i < j],
    node_color="red",
    edge_color="blue",
)
print(len(sol))
15

クリーク被覆問題

無向グラフを,最小の数のクリークで被覆する問題がクリーク被覆問題(clique cover problem)である. 点を被覆する問題(クリーク点被覆問題)は,補グラフを作ればグラフ彩色問題に帰着できる.枝の被覆問題(クリーク枝被覆問題)もあり, 両者とも $NP$-困難である.

極大クリークを列挙すると,両者ともに集合被覆問題に帰着できる.以下に,点被覆の例を示す.

model = Model()
x = {}
nodes = {}
for i, c in enumerate(nx.find_cliques(G)):
    x[i] = model.addVar(vtype="B", name=f"x[{i}]")
    nodes[i] = c
cliques = {}  # 点iを含んでいるクリークの集合
for i in range(n):
    cliques[i] = []
for c in nodes:
    for i in nodes[c]:
        cliques[i].append(c)
for i in range(n):
    model.addConstr(quicksum(x[c] for c in cliques[i]) >= 1)
model.setObjective(quicksum(x[c] for c in x), GRB.MINIMIZE)
Academic license - for non-commercial use only - expires 2022-04-03
Using license file /Users/mikiokubo/gurobi.lic
model.optimize()
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 1000 rows, 55758 columns and 4183991 nonzeros
Model fingerprint: 0x4c062606
Variable types: 0 continuous, 55758 integer (55758 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 30.0000000
Presolve removed 17 rows and 364 columns (presolve time = 5s) ...
Presolve removed 68 rows and 364 columns (presolve time = 10s) ...
Presolve removed 70 rows and 2762 columns (presolve time = 16s) ...
Presolve removed 107 rows and 2763 columns (presolve time = 20s) ...
Presolve removed 107 rows and 2763 columns (presolve time = 25s) ...
Presolve removed 107 rows and 2763 columns (presolve time = 30s) ...
Presolve removed 107 rows and 2763 columns
Presolve time: 30.61s
Presolved: 893 rows, 52995 columns, 3958930 nonzeros
Variable types: 0 continuous, 52995 integer (52995 binary)

Root simplex log...

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   8.930000e+02   0.000000e+00     31s

Starting sifting (using dual simplex for sub-problems)...

    Iter     Pivots    Primal Obj      Dual Obj        Time
       0          0     infinity      0.0000000e+00     31s
       1        171   2.3000000e+01   5.3076923e+00     31s
       2        347   2.1000000e+01   7.4000000e+00     31s
       3        821   1.9500000e+01   9.3033708e+00     31s
       4       1681   1.9000000e+01   1.1172336e+01     31s
       5       3300   1.7772358e+01   1.2865058e+01     32s
       6       6533   1.6201641e+01   1.3595413e+01     32s
       7       8637   1.6057971e+01   1.4183650e+01     33s
       8      11303   1.6000000e+01   1.4409530e+01     33s
       9      14098   1.6000000e+01   1.4582271e+01     34s
      10      15686   1.6000000e+01   1.4724554e+01     34s
      11      17102   1.6000000e+01   1.4881946e+01     35s
      12      18681   1.6000000e+01   1.4971198e+01     35s

Sifting complete

   21452    1.6000000e+01   0.000000e+00   0.000000e+00     36s

Root relaxation: objective 1.600000e+01, 21452 iterations, 5.33 seconds
Total elapsed time = 40.43s

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   16.00000    0  123   30.00000   16.00000  46.7%     -   41s
H    0     0                      18.0000000   16.00000  11.1%     -   41s
     0     0   16.00000    0  167   18.00000   16.00000  11.1%     -   46s
H    0     0                      17.0000000   16.00000  5.88%     -   52s
     0     0   16.00000    0  138   17.00000   16.00000  5.88%     -   58s
     0     0   16.00000    0  120   17.00000   16.00000  5.88%     -   83s
     0     0   16.00000    0  169   17.00000   16.00000  5.88%     -   91s
     0     0   16.00000    0  142   17.00000   16.00000  5.88%     -   93s
     0     0   16.00000    0  128   17.00000   16.00000  5.88%     -   99s
     0     0   16.00000    0   98   17.00000   16.00000  5.88%     -  105s
     0     0   16.00000    0  117   17.00000   16.00000  5.88%     -  111s
     0     0   16.00000    0  135   17.00000   16.00000  5.88%     -  113s
     0     0   16.00000    0  119   17.00000   16.00000  5.88%     -  118s
     0     0   16.00000    0  116   17.00000   16.00000  5.88%     -  125s
     0     0   16.00000    0  116   17.00000   16.00000  5.88%     -  127s
     0     0   16.00000    0  116   17.00000   16.00000  5.88%     -  129s
     0     2   16.00000    0  116   17.00000   16.00000  5.88%     -  132s
     1     3   16.00000    1  176   17.00000   16.00000  5.88%  3306  135s
     7     2   16.00000    4  168   17.00000   16.00000  5.88%   846  140s

Cutting planes:
  MIR: 4

Explored 11 nodes (69694 simplex iterations) in 141.09 seconds
Thread count was 16 (of 16 available processors)

Solution count 3: 17 18 30 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.700000000000e+01, best bound 1.700000000000e+01, gap 0.0000%
lhs = np.zeros(n)
for c in x:
    if x[c].X > 0.001:
        for i in nodes[c]:
            lhs[i] += 1
edges = []
for c in x:
    if x[c].X > 0.001:
        # print(nodes[c])
        for i in nodes[c]:
            for j in nodes[c]:
                if i != j:
                    edges.append((i, j))
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=10, edge_color="blue")
nx.draw_networkx_edges(G, pos=pos, edgelist=edges, edge_color="orange", width=2)
plt.show()