グラフの埋め込み問題をQUBOに変換してGurobiで解く

# ポイント
* ゲストグラフ$G$をホストグラフ$H$に埋め込むグラフの埋め込み問題を等価なQUBOに変換する方法
* PyQUBOによるQUBOの生成
* GurobiによるQUBOの解を求める方法
* 埋め込みが成功したら，解探索を終了するGurobiのコールバック関数の設定法
* NetworkXを用いた無向グラフの描画

# グラフの埋め込み問題
問題の入力は，ゲストグラフ$G=(V_G,E_G)$とホストグラフ$H=(V_H,E_H)$です．いずれも無向グラフで$V_G$と$V_H$が頂点集合，$E_G$と$E_H$が辺集合，つまり頂点の無順序対です．求めるのは$V_G$から$V_H$への**埋め込み**，つまり写像$f:V_G\rightarrow V_H$です．ゲストグラフの辺集合$E_G$に写像$f$を適用して得られる辺集合$E_{f(G)}=\lbrace (f(u),f(v))| (u,v)\in E_G\rbrace$がすべてホストグラフの辺集合$E_H$に含まれている，つまり$E_{f(G)}\subseteq E_H$であれば，埋め込み成功です．このとき，グラフ$H$は$GG$を部分グラフとして含みます．埋め込みできない場合は，埋め込みができた辺集合$E_{f(G)}\cap E_H$の要素数を最大化するように埋め込み$f$を求めます．

## グラフ埋め込みの解の例
下記の図では，8頂点からなるゲストグラフ$G$を，9頂点のホストグラフ$H$に埋め込む例を示しています．$G$は12辺で構成されており，埋め込みに成功した辺は青色で，失敗した辺は赤色で示されています．$H$では，埋め込みに利用された辺が青色で表示されています．12辺のうち，埋め込みに成功したのは9辺でした．

|ゲストグラフ$G$: $2\times 2 \times 2$のキューブ|ホストグラフ$H$: $3\times 3$のグリッド|
|---|---|
|<img width="300" alt="guest2x2x2.png" src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3742233/ae9f457a-75dd-e218-59b8-0009510a20d4.png>|<img width="300" alt="host3x3.png" src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3742233/eefcefb8-af3b-1e1f-278c-740c85b4e8c0.png>|

# グラフ埋め込み問題のQUBOへの変換
ゲストグラフ$G=(V_G,E_G)$とホストグラフ$H=(V_H,E_H)$において，QUBOのバイナリ変数の集合$X$を以下のように定義します．
$$
X=\lbrace x_{u,v} | u\in V_G, v\in V_H\rbrace
$$
ここで，埋め込みを示す写像$f:V_G\rightarrow V_H$が$f(u)=v$を満たす場合，$x_{u,v}=1$となるようにQUBOを設計します．
まず，任意の$u\in V_G$に対して，$f(u)$が一意に定まるためように，以下の式を用います．
$$
\sum_{v\in V_H} x_{u,v} =1
$$
さらに，すべての$v\in V_H$について，$f(u)=v$となる$u$が存在しないか，1つだけ存在するため，次の式を満たす必要があります．
$$
\sum_{u\in V_G} x_{u,v} = 0 \quad または \quad 1
$$
前者はone-hot制約[^1]，後者はzero-one-hot制約[^1]を満たすように，QUBOの式を以下のように設定します．
$$
C(X)=\sum_{u\in V_G}(1-\sum_{v\in V_H} x_{u,v})^2+\sum_{v\in V_H}(1-\sum_{u\in V_G} x_{u,v})\sum_{u\in V_G}x_{u,v}
$$
$C(X)$が最小値0を取る場合，$X=(x_{u,v})$は上記の制約をそれぞれ満たし，$x_{u,v}=1$は$f(u)=v$を示します．
次に、成功した埋め込みの辺の総数を計算するQUBOを構築します．$E_G$と$E_H$からそれぞれ一つずつ辺を取り，$(u,u')\in E_G$，$(v,v')\in E_H$とします．埋め込み$f$が$f(u)=v$かつ$f(u')=v'$，または$f(u)=v'$かつ$f(u')=v$を満たす場合，辺$(u,u')$は$(v,v')$に埋め込まれています．したがって，埋め込まれた辺の総数は次の式で算出されます．
$$
E(X)=\sum_{(u,u')\in E_G}\sum_{(v,v')\in E_H} (x_{u,v}x_{u',v'}+x_{u,v'}x_{u',v})
$$
最終的に，$C(X)$が最小値0を取り，かつ$E(X)$が最大値を取るようにQUBOを設定します．
$$
  P\cdot C(X)-E(X)
$$
ここで，$P$は$C(X)$を$E(X)$より優先するための係数で，ゲストグラフ$G$の頂点数，つまり$P=|V_G|$とすれば十分です．

# Pythonプログラム
ホストグラフ$G$とゲストグラフ$H$は、NetworkXの無向グラフとしてそれぞれ`host`と`guest`に指定されます．この例では，`guest`には$2\times 2 \times 2$のキューブ，`host`には$3\times 3$のグリッドが設定されていますが，これらを変更することで，異なるグラフの埋め込み問題を解くことが可能です.

関数`embedding_qubo`は、PyQUBOを使用して$P\cdot C(X)-E(X)$の式を`formula`として構築し、これを展開してQUBO形式の`qubo`に変換します。この`qubo`には定数項が含まれず、定数項は`offset`として別途保持されます。コード中では、`guest_1+host_01`が$C(X)$を表し、`match`は$-E(X)$に相当します。

関数`gurobi`は、`qubo`の解`solution`を求めます．`vars`は、キーとしてタプル$(u,v)$の文字列を持ち、値として$x_{u,v}$の値を持つ辞書です．`qubo`をGurobiの目的関数に適合する形式で`obj`に変換しています．定数項`offset`とゲストグラフの辺の数`len(guest.edges)`の合計に対応するマイナスの値を持つ解が見つかった場合，解探索を即座に終了させるためのコールバック関数`mycallback`を設定しています．この条件が満たされた場合，すべてのゲストグラフの辺が埋め込みに成功していることになります．

辞書`mapping`は，キーとして$u$，値として$f(u)$を持ち，`solution`から導出されます.
$x_{u,v}=1$のときに、キー$u$に対する値として$v$が割り当てられます．

関数`print_embedding`は，`solution`をわかりやすくテキスト表示するためのもので、`draw_embedding`関数はmatplotlibを利用して、`guest`から`host`への埋め込み結果を視覚的に描画し，画像ファイルに保存します．


```Python
#!/usr/bin/env python3

import argparse
import networkx as nx
import matplotlib.pyplot as plt
import pyqubo as pq
import gurobipy as gp
import ast


# generate a QUBO model for embedding guest graph in host graph
def embedding_qubo(guest, host):
    pairs = [(v_g, v_h) for v_h in host.nodes for v_g in guest.nodes]
    vars = {p: pq.Binary(str(p)) for p in pairs}
    host2guest = [[vars[(v_g, v_h)] for v_g in guest.nodes] for v_h in host.nodes]
    guest2host = [[vars[(v_g, v_h)] for v_h in host.nodes] for v_g in guest.nodes]
    host_01 = sum([sum(p) * (sum(p) - 1) for p in host2guest])
    guest_1 = sum([(sum(p) - 1) ** 2 for p in guest2host])
    match = sum(
        [
            -vars[(e_g[0], e_h[0])] * vars[(e_g[1], e_h[1])]
            - vars[(e_g[0], e_h[1])] * vars[(e_g[1], e_h[0])]
            for e_h in host.edges
            for e_g in guest.edges
        ]
    )
    formula = match + len(guest.nodes) * host_01 + len(guest.nodes) * guest_1
    model = formula.compile()
    qubo, offset = model.to_qubo()
    return qubo, int(offset)


# Gurobi call back function: terminate if all guest edges are embedded
def mycallback(model, where):
    if where == gp.GRB.Callback.MIPSOL:
        if int(model.cbGet(gp.GRB.Callback.MIPSOL_OBJ)) == qubo_optimal:
            model.terminate()


# Find solution of qubo using Gurobi
def gurobi(qubo, timelimit):
    model = gp.Model()
    model.setParam("TimeLimit", timelimit)
    vars = {}
    for v in list((v for tup in qubo for v in tup)):
        vars[v] = model.addVar(vtype=gp.GRB.BINARY, name=v)
    obj = [val * vars[key[0]] * vars[key[1]] for key, val in qubo.items()]
    model.setObjective(gp.quicksum(obj), sense=gp.GRB.MINIMIZE)
    model.optimize(mycallback)
    solution = {key: int(val.X) for key, val in vars.items()}
    return solution, int(model.ObjVal)


# Display embedding result
def print_embedding(guest, host, mapping):
    print("Node mapping: guest to host")
    for key, val in sorted(mapping.items()):
        print(key, ":", val)
    print("Edge mapping: guest to host")
    for u, v in guest.edges:
        if (mapping[u], mapping[v]) in host.edges:
            print((u, v), ":", (mapping[u], mapping[v]))
        else:
            print((u, v), ": failure")


# Draw embedding
def draw_embedding(guest, host, mapping, guest_file, host_file):
    combined = nx.Graph()
    combined.add_edges_from([(mapping[u], mapping[v]) for u, v in guest.edges])
    edge_color = ["blue" if e_h in combined.edges else "black" for e_h in host.edges]
    node_color = [
        "lightblue" if v_h in mapping.values() else "gray" for v_h in host.nodes
    ]

    pos = nx.spring_layout(host)
    nx.draw(
        host,
        pos,
        with_labels=True,
        node_color=node_color,
        edge_color=edge_color,
        width=2,
    )
    plt.savefig(host_file)
    plt.clf()
    edge_color = [
        "blue" if (mapping[u], mapping[v]) in host.edges else "red"
        for u, v in guest.edges
    ]
    pos = nx.spring_layout(host)
    plt.clf()
    nx.draw(
        guest,
        with_labels=True,
        node_color="lightblue",
        edge_color=edge_color,
        width=2,
    )
    plt.savefig(guest_file)


def main():
    global qubo_optimal
    parser = argparse.ArgumentParser(
        description="Embedding a guest graph in a host graph using Pyqubo and Gurobi"
    )
    parser.add_argument(
        "-t", "--timelimit", default=30, type=int, help="Gurobi Timelimit"
    )
    parser.add_argument(
        "--guest", default="guest.png", type=str, help="File for drawing guest graph"
    )
    parser.add_argument(
        "--host", default="host.png", type=str, help="File for drawing host graph"
    )

    args = parser.parse_args()
    host = nx.grid_graph([3, 3])  # host graph
    guest = nx.grid_graph([2, 2, 2])  # guest graph to be embedded

    qubo, offset = embedding_qubo(guest, host)  # generate QUBO for embedding
    qubo_optimal = -len(guest.edges) - offset  # optimal QUBO energy if all guest edge
    solution, objval = gurobi(qubo, args.timelimit)  # Find solution using Gurobi
    mapping = {  # Find guest to host node mapping from solution
        ast.literal_eval(key)[0]: ast.literal_eval(key)[1]
        for key, val in solution.items()
        if val == 1
    }
    print_embedding(guest, host, mapping)  # display obtained embedding
    draw_embedding(guest, host, mapping, args.guest, args.host)  # draw embedding
    print(f"optimal = {-len(guest.edges)} obtained = {offset+objval}")


if __name__ == "__main__":
    main()
```

# 実行結果
先のPythonプログラムの実行結果です．ゲストグラフ$G$のノードは，$(i,j,k)$ ($0\leq i,j,k\leq 1$)で表され，ホストグラフ$H$のノードは$(i,j)$ ($0\leq i,j\leq 2$)で表されています．ゲストグラフ$G$のうち3つの辺の埋め込みが失敗しています．Gurobiはこの解がBoundと一致する最適解で，より良い埋め込みは存在しないことを示しています．

```
Node mapping: guest to host
(0, 0, 0) : (2, 2)
(0, 0, 1) : (1, 0)
(0, 1, 0) : (2, 1)
(0, 1, 1) : (1, 1)
(1, 0, 0) : (1, 2)
(1, 0, 1) : (0, 0)
(1, 1, 0) : (0, 2)
(1, 1, 1) : (0, 1)
Edge mapping: guest to host
((0, 0, 0), (1, 0, 0)) : ((2, 2), (1, 2))
((0, 0, 0), (0, 1, 0)) : ((2, 2), (2, 1))
((0, 0, 0), (0, 0, 1)) : failure
((0, 0, 1), (1, 0, 1)) : ((1, 0), (0, 0))
((0, 0, 1), (0, 1, 1)) : ((1, 0), (1, 1))
((0, 1, 0), (1, 1, 0)) : failure
((0, 1, 0), (0, 1, 1)) : ((2, 1), (1, 1))
((0, 1, 1), (1, 1, 1)) : ((1, 1), (0, 1))
((1, 0, 0), (1, 1, 0)) : ((1, 2), (0, 2))
((1, 0, 0), (1, 0, 1)) : failure
((1, 0, 1), (1, 1, 1)) : ((0, 0), (0, 1))
((1, 1, 0), (1, 1, 1)) : ((0, 2), (0, 1))
optimal = -12 obtained = -9
```


[^1]:[N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解く](https://qiita.com/koji_nakano/items/860bd9fb46072af986af)

# 参考
* [Various QUBO Formulations of the Graph Isomorphism Problem and Related Problems](http://www.bncss.org/index.php/bncss/article/view/178)