最大クリーク問題の定式化

# これなに
2017/10/7 の[第3回ORセミナー](http://www.orsj.or.jp/activity/seminar.html#semi3)で「最大クリーク問題の定式化」の質問がありましたので、やってみました。

## 最大安定集合問題と最大クリーク問題の関係

- [最大安定問題](https://qiita.com/SaitoTsutomu/items/33ad03bcaa421eb2ba17)は、選ばれた点間で**辺がない**。
- 最大クリーク問題は、選ばれた点間で**辺がある**。

このことから、最大安定問題の最適解と「その[補グラフ](https://ja.wikipedia.org/wiki/%E8%A3%9C%E3%82%B0%E3%83%A9%E3%83%95)[^1]の最大クリーク問題」の最適解は同じものにすることができます。

[^1]: 補グラフとは、辺の有無を逆転させたグラフ。

## 最大安定集合問題の定式化

セミナーでも紹介されていましたが、定式化は、「任意の辺の両端のうち1つまで選べる」という制約条件だけでシンプルです。

最大化|$\sum_{i \in 点集合}{~~~~ x_i}$ (1)
:--|:--
制約条件|$x_i + x_j \le 1 ~~~ \forall i,j \in 辺集合$ (2)
 |$x_i \in \\{0,1\\}~~~ \forall i \in 点集合$ (3)

### pythonによる実装

```py3:python
from ortoolpy import maximum_stable_set
maximum_stable_set??
>>>
def maximum_stable_set(g, weight='weight'):
    """
    最大安定集合問題
    入力
        g: グラフ(node:weight)
        weight: 重みの属性文字
    出力
        最大安定集合の重みの合計と頂点番号リスト
    """
    from pulp import LpProblem, LpMaximize, LpBinary, lpDot, lpSum, value
    m = LpProblem(sense=LpMaximize)
    v = [addvar(cat=LpBinary) for _ in g.nodes()] # (3)
    for i, j in g.edges():
        m += v[i] + v[j] <= 1 # (2)
    m += lpDot([g.node[i].get(weight, 1) for i in g.nodes()], v) # (1)
    if m.solve() != 1: return None
    return value(m.objective), [i for i, x in enumerate(v) if value(x) > 0.5]
```

## 最大クリーク問題の定式化

選んだ点に応じて、クリーク(完全グラフ)になるように両端が選ばれた辺の数を強制することにより、定式化できます(実際に解く場合は、最大安定集合問題を解いた方がよいです)。

最大化|$S_v$ (1)
:--|:--
制約条件|$E_{ij} \le V_i ~~~ \forall i,j \in 辺集合$ (2)
 |$E_{ij} \le V_j ~~~ \forall i,j \in 辺集合$ (2)
 |$S_v = \sum_{i \in 点集合}{~~~~ V_i}$ (3)
 |$S_e = \sum_{i,j \in 辺集合}{~~~~ E_{ij}}$ (4)
 |$S_e$は、折れ線 curve の $S_v$ の位置以上 (5)
 |$V_i \in \\{0,1\\}~~~ \forall i \in 点集合$ (6)
 |$E_{ij} \in \\{0,1\\}~~~ \forall i,j \in 辺集合$ (7)
 |$S_v, S_e \ge 0$ (8)


### pythonによる実装

```py3:python
import networkx as nx
from itertools import accumulate, count, takewhile
from pulp import *
from ortoolpy import addlines_conv
def maximum_clique(g):
    """
    最大クリーク問題
    入力
        g: グラフ
    出力
        最大クリークの頂点番号リスト
    """
    m = LpProblem(sense=LpMaximize) # 数理モデル
    for i in g.nodes():
        g.node[i]['VarV'] = LpVariable(f'V{i}', cat=LpBinary) # 点変数(6)
    for i,j in g.edges():
        g[i][j]['VarE'] = LpVariable(f'E{i}_{j}', cat=LpBinary) # 辺変数(7)
        m += g[i][j]['VarE'] <= g.node[i]['VarV'] # (2)
        m += g[i][j]['VarE'] <= g.node[j]['VarV'] # (2)
    vsv = LpVariable('VarSV') # 点数変数 # (8)
    vse = LpVariable('VarSE') # 辺数変数 # (8)
    m += vsv == lpSum(g.node[i]['VarV'] for i in g.nodes()) # (3)
    m += vse == lpSum(g[i][j]['VarE'] for i,j in g.edges()) # (4)
    m += vsv # 目的関数(1)
    curve = list(takewhile(lambda x:x[1]<=g.number_of_edges(),
        enumerate(accumulate(count()),1)))
    addlines_conv(m, curve, vsv, vse) # 選んだサブグラフをクリークに強制(5)
    m.solve()
    return [i for i in g.nodes() if value(g.node[i]['VarV'])>0.5]
```

### 補足
ここで curveは完全グラフの点数と辺数に対応します。
`addlines_conv(m, curve, vsv, vse)`は「vse >= 『curveのvsv(点数)に対応する値(辺数)』」となる制約条件を数理モデルmに追加します。
例えば、選択した点の数であるvsvが4だと、vseは6以上にしなければいけません。vseは$E_{ij}$の和です。また、$E_{ij}$は制約条件の(2)から、両端点を選択しないと1にできません。すなわち4点かつ6辺以上のグラフは、辺数が6のクリークになります。

```py3:python
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
curve = list(takewhile(lambda x:x[1]<=28,
    enumerate(accumulate(count()),1)))
print(curve)
plt.plot(*np.array(curve).T);
>>>
[(1, 0), (2, 1), (3, 3), (4, 6), (5, 10), (6, 15), (7, 21), (8, 28)]
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/13955/ebeabdc6-08b3-98fd-3fec-4f47084bab51.png)

### サンプルで確認

最大クリークが得られるか確認します。

```py3:python
g = nx.fast_gnp_random_graph(6, 0.3, 3)
nx.draw_networkx(g)
maximum_clique(g)
>>>
[2, 3, 5]
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/13955/8006ecba-7db4-086f-d5c0-0b85e3146bb1.png)

もう1つ確認。

```py3:python
g = nx.fast_gnp_random_graph(7, 0.4, 49)
nx.draw_networkx(g)
maximum_clique(g)
>>>
[0, 3, 4, 5]
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/13955/39e8c197-9d24-85dd-8e87-f49660fa09b2.png)

### 追記
「点を選んだら、vsv本の辺と接続すること」という定式化の方が自然かもしれません。

以上