# 部分グラフ同型問題 :::{container} prog-cpp 2つの無向グラフ $G_H=(V_H,E_H)$（ホストグラフ）と $G_G=(V_G,E_G)$（ゲストグラフ）が与えられたとき,**部分グラフ同型問題**は $G_H$ が $G_G$ と同型な部分グラフを含むかどうかを判定する問題です. より形式的には,すべての辺 $(u,v)\in E_G$ に対して $(\sigma(u),\sigma(v))$ がホストグラフの辺でもある（すなわち $(\sigma(u),\sigma(v))\in E_H$）ような**単射** $\sigma:V_G\rightarrow V_H$ を見つけることが目標です. 例として,以下のホストグラフとゲストグラフを考えます: 
10頂点のホストグラフ $G_H=(V_H,E_H)$ の例

6頂点のゲストグラフ $G_G=(V_G,E_G)$ の例

解 $\sigma$ の一例は次の通りです: | $G_G$ の頂点 $i$ | 0 | 1 | 2 | 3 | 4 | 5 | |:---:|:---:|:---:|:---:|:---:|:---:|:---:| | $G_H$ の頂点 $\sigma(i)$ | 1 | 4 | 6 | 7 | 9 | 8 | この解は次のように可視化されます: 
部分グラフ同型問題の解

## 部分グラフ同型問題のQUBO定式化 **ゲストグラフ** $G_G=(V_G,E_G)$ が $m$ 個の頂点（ラベル $0, 1, \ldots m-1$）を持ち,**ホストグラフ** $G_H=(V_H,E_H)$ が $n$ 個の頂点（ラベル $0, 1, \ldots n-1$）を持つとします. $mn$ 個のバイナリ変数を持つ $m\times n$ の**バイナリ行列** $X=(x_{i,j})$（$0\leq i\leq m-1, 0\leq j\leq n-1$）を導入します. この行列は単射 $\sigma:V_G\rightarrow V_H$ を表し,$x_{i,j}=1$ は $\sigma(i)=j$ の場合です. 例えば,部分グラフ同型問題の解は以下の $6\times 10$ バイナリ行列で表現できます: | $i$ | $\sigma(i)$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:| | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | | 1 | 4 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | | 2 | 6 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | | 3 | 7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 4 | 9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 5 | 8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | $X$ は単射を表すため,以下の制約を満たす必要があります: - **行制約**: 各ゲスト頂点はちょうど1つのホスト頂点に写像される.すなわち各行の和が1. - **列制約**: 各ホスト頂点は高々1つのゲスト頂点から写像される.すなわち各列の和が0または1. これらをまとめると,すべての制約が満たされたときに最小値をとる以下の**Hi-QUBO形式の制約**になります: $$ \begin{aligned} \text{constraint} &= \sum_{i=0}^{m-1}\Bigl(\sum_{j=0}^{n-1}x_{i,j} = 1\Bigr)+\sum_{j=0}^{n-1}\Bigl(0\leq \sum_{i=0}^{m-1}x_{i,j} \leq 1\Bigr) \end{aligned} $$ QUBO形式では,同じ制約を次のように表現できます: $$ \begin{aligned} \text{constraint} &= \sum_{i=0}^{m-1}\Bigl(\sum_{j=0}^{n-1}x_{i,j} - 1\Bigr)^2+\sum_{j=0}^{n-1}\sum_{i=0}^{m-1}x_{i,j}\Bigl(\sum_{i=0}^{m-1}x_{i,j}-1\Bigr) \end{aligned} $$ 次に,目的関数をホスト辺に写像されたゲスト辺の数として定義します: $$ \begin{aligned} \text{objective} &= \sum_{(u_G,v_G)\in E_G}\sum_{(u_H,v_H)\in E_H} (x_{u_G,u_H}x_{v_G,v_H}+x_{u_G,v_H}x_{v_G,u_H}) \end{aligned} $$ ここで,無向のゲスト辺 $(u_G,v_G)\in E_G$ はホスト辺 $(u_H,v_H)\in E_H$ に2つの対称的な方法で対応できます: - $(u_G, v_G)\mapsto (u_H,v_H)$ - $(u_G, v_G)\mapsto (v_H,u_H)$ したがって,2次の項 $x_{u_G,u_H}x_{v_G,v_H}$ と $x_{u_G,v_H}x_{v_G,u_H}$ の両方を含めます. 最終的に,目的関数と制約を1つのQUBO式にまとめます: $$ \begin{aligned} f &= -\text{objective} + mn\times \text{constraint} \end{aligned} $$ ペナルティ係数 $mn$ は,目的関数の改善よりも制約の充足を優先するために選ばれています. $f$ の最良値は,制約項が0で目的関数がゲスト辺の数に等しいときに達成されます. ## 部分グラフ同型問題のHi-QUBOプログラム 上記のQUBO定式化に基づき,以下のHi-QUBOプログラムは $M=6$ 頂点のゲストグラフと $N=10$ 頂点のホストグラフに対する部分グラフ同型問題を解きます: ```{literalinclude} /../programFiles/cppPrograms/example/graph/subgraph-isomorphism-program1.cpp :language: cpp :caption: subgraph-isomorphism-program1.cpp ``` ゲストグラフとホストグラフは,それぞれ辺リスト `guest` と `host` として与えられます. $M\times N$ のバイナリ行列 `x` を定義し,上記の定式化に従って `constraint`,`objective`,`f` を構成します. Easy Solver のインスタンスを `f` に対して作成し,目標エネルギーを $-|E_G|$（ゲスト辺数の負の値）に設定します.これはすべてのゲスト辺がホスト辺に写像されたときの `-objective` の最良値です. 得られた解は `sol` に格納されます. `sol` の下での `x`,`objective`,`constraint` の値が出力されます. 関数 **`qbpp::onehot_to_int()`** を用いて,ゲスト頂点からホスト頂点への写像（`guest_to_host`,$\sigma$）とホスト頂点からゲスト頂点への写像（`host_to_guest`,$\sigma^{-1}$）も出力します. ゲストグラフとホストグラフはそれぞれ `guest_graph.svg` と `host_graph.svg` として保存されます. 最後に,解が `subgraph_isomorphism.svg` に可視化されます.写像で選択されたホスト頂点と,ゲスト辺に対応するホスト辺がハイライトされています. このプログラムは以下の出力を生成します: ```{include} /../programFiles/markDown/example/graph/subgraph-isomorphism.md :start-after: :end-before: ``` 目的関数値はゲスト辺の数（$|E_G|=8$）に等しく,すべての制約が満たされています（`constraint` = 0）. したがって,プログラムは有効な部分グラフ同型に対応する最適解を見つけました. host_to_guest のエントリが `-1` の場合,対応するホスト頂点にはゲスト頂点が写像されていないことを意味します. ## `einsum` を使った簡潔な目的関数 `objective` を構築する二重の辺ループは,各グラフを辺リストではなく **二値の隣接行列**として表現することで,[`qbpp::einsum`](../../advanced/Einstein-sum.md) を使って 書き換えられます.次のように上三角の隣接行列を用意します: - `A_G[a, b] = 1` ($\{a, b\} \in E_G$ のとき,無向ゲスト辺ごとに 1 エントリ) - `A_H[c, d] = 1` ($\{c, d\} \in E_H$ のとき,無向ホスト辺ごとに 1 エントリ) すると目的関数 $\sum_{(u_G,v_G)\in E_G}\sum_{(u_H,v_H)\in E_H} (x_{u_G,u_H}x_{v_G,v_H}+x_{u_G,v_H}x_{v_G,u_H})$ は 2 つの `einsum` 呼び出しの和として書けます: ```{include} /../programFiles/markDown/example/graph/subgraph-isomorphism.md :start-after: :end-before: ``` subscript `"ab,cd,ac,bd->"` はそのまま $\sum_{a,b,c,d} A_G[a,b]\, A_H[c,d]\, x_{a,c}\, x_{b,d}$ を表しており, [`einsum` ドキュメント](../../advanced/Einstein-sum.md)の QAP 形のテンソル縮約と同じパターンです. 2 つ目の呼び出しはホスト軸を入れ替えた $(u_G, v_G) \mapsto (v_H, u_H)$ の対称写像をカバーします （`ac,bd` の代わりに `ad,bc`）. 得られる QUBO 式の整理後の項集合は for ループ版と完全に同じですが, 構築コードは大幅に短くなり,`einsum` 内部でマルチスレッドによる並列化も 効きます.トレードオフはメモリ使用量で,辺リスト方式は $|E_G|+|E_H|$ に 比例するのに対し,隣接行列表現は $\Theta(M^2 + N^2)$ となるため, 非常に疎で巨大なグラフでは for ループ版の方が依然として有利です. ::: :::{container} prog-python 2つの無向グラフ $G_H=(V_H,E_H)$（ホストグラフ）と $G_G=(V_G,E_G)$（ゲストグラフ）が与えられたとき,**部分グラフ同型問題**は $G_H$ が $G_G$ と同型な部分グラフを含むかどうかを判定する問題です. より形式的には,すべての辺 $(u,v)\in E_G$ に対して $(\sigma(u),\sigma(v))$ がホストグラフの辺でもある（すなわち $(\sigma(u),\sigma(v))\in E_H$）ような**単射** $\sigma:V_G\rightarrow V_H$ を見つけることが目標です. 例として,以下のホストグラフとゲストグラフを考えます: 
10頂点のホストグラフ $G_H=(V_H,E_H)$ の例

6頂点のゲストグラフ $G_G=(V_G,E_G)$ の例

解 $\sigma$ の一例は次の通りです: | $G_G$ の頂点 $i$ | 0 | 1 | 2 | 3 | 4 | 5 | |:---:|:---:|:---:|:---:|:---:|:---:|:---:| | $G_H$ の頂点 $\sigma(i)$ | 1 | 4 | 6 | 7 | 9 | 8 | この解は次のように可視化されます: 
部分グラフ同型問題の解

## 部分グラフ同型問題のQUBO定式化 **ゲストグラフ** $G_G=(V_G,E_G)$ が $m$ 個の頂点（ラベル $0, 1, \ldots m-1$）を持ち,**ホストグラフ** $G_H=(V_H,E_H)$ が $n$ 個の頂点（ラベル $0, 1, \ldots n-1$）を持つとします. $mn$ 個のバイナリ変数を持つ $m\times n$ の**バイナリ行列** $X=(x_{i,j})$（$0\leq i\leq m-1, 0\leq j\leq n-1$）を導入します. この行列は単射 $\sigma:V_G\rightarrow V_H$ を表し,$x_{i,j}=1$ は $\sigma(i)=j$ の場合です. 例えば,部分グラフ同型問題の解は以下の $6\times 10$ バイナリ行列で表現できます: | $i$ | $\sigma(i)$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | |:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:| | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | | 1 | 4 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | | 2 | 6 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | | 3 | 7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | | 4 | 9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 5 | 8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | $X$ は単射を表すため,以下の制約を満たす必要があります: - **行制約**: 各ゲスト頂点はちょうど1つのホスト頂点に写像される.すなわち各行の和が1. - **列制約**: 各ホスト頂点は高々1つのゲスト頂点から写像される.すなわち各列の和が0または1. これらをまとめると,すべての制約が満たされたときに最小値をとる以下の**Hi-QUBO形式の制約**になります: $$ \begin{aligned} \text{constraint} &= \sum_{i=0}^{m-1}\Bigl(\sum_{j=0}^{n-1}x_{i,j} = 1\Bigr)+\sum_{j=0}^{n-1}\Bigl(0\leq \sum_{i=0}^{m-1}x_{i,j} \leq 1\Bigr) \end{aligned} $$ QUBO形式では,同じ制約を次のように表現できます: $$ \begin{aligned} \text{constraint} &= \sum_{i=0}^{m-1}\Bigl(\sum_{j=0}^{n-1}x_{i,j} - 1\Bigr)^2+\sum_{j=0}^{n-1}\sum_{i=0}^{m-1}x_{i,j}\Bigl(\sum_{i=0}^{m-1}x_{i,j}-1\Bigr) \end{aligned} $$ 次に,目的関数をホスト辺に写像されたゲスト辺の数として定義します: $$ \begin{aligned} \text{objective} &= \sum_{(u_G,v_G)\in E_G}\sum_{(u_H,v_H)\in E_H} (x_{u_G,u_H}x_{v_G,v_H}+x_{u_G,v_H}x_{v_G,u_H}) \end{aligned} $$ ここで,無向のゲスト辺 $(u_G,v_G)\in E_G$ はホスト辺 $(u_H,v_H)\in E_H$ に2つの対称的な方法で対応できます: - $(u_G, v_G)\mapsto (u_H,v_H)$ - $(u_G, v_G)\mapsto (v_H,u_H)$ したがって,2次の項 $x_{u_G,u_H}x_{v_G,v_H}$ と $x_{u_G,v_H}x_{v_G,u_H}$ の両方を含めます. 最終的に,目的関数と制約を1つのQUBO式にまとめます: $$ \begin{aligned} f &= -\text{objective} + mn\times \text{constraint} \end{aligned} $$ ペナルティ係数 $mn$ は,目的関数の改善よりも制約の充足を優先するために選ばれています. $f$ の最良値は,制約項が0で目的関数がゲスト辺の数に等しいときに達成されます. ## 部分グラフ同型問題のPyQBPPプログラム 上記のQUBO定式化に基づき,以下のPyQBPPプログラムは $M=6$ 頂点のゲストグラフと $N=10$ 頂点のホストグラフに対する部分グラフ同型問題を解きます: ```{literalinclude} /../programFiles/pythonPrograms/example/graph/subgraph-isomorphism-program1.py :language: python :caption: subgraph-isomorphism-program1.py ``` ゲストグラフとホストグラフは,それぞれ辺リスト `guest` と `host` として与えられます. $M\times N$ のバイナリ行列 `x` を定義し,上記の定式化に従って `constraint`,`objective`,`f` を構成します. Easy Solver のインスタンスを `f` に対して作成し,目標エネルギーを $-|E_G|$（ゲスト辺数の負の値）に設定します.これはすべてのゲスト辺がホスト辺に写像されたときの `-objective` の最良値です. 得られた解は `sol` に格納されます. `sol` の下での `x`,`objective`,`constraint` の値が出力されます. 関数 **`qbpp.onehot_to_int()`** を用いて,ゲスト頂点からホスト頂点への写像（`guest_to_host`,$\sigma$）とホスト頂点からゲスト頂点への写像（`host_to_guest`,$\sigma^{-1}$）も出力します. このプログラムは以下の出力を生成します: ```{include} /../programFiles/markDown/example/graph/subgraph-isomorphism.md :start-after: :end-before: ``` 目的関数値はゲスト辺の数（$|E_G|=8$）に等しく,すべての制約が満たされています（`constraint` = 0）. したがって,プログラムは有効な部分グラフ同型に対応する最適解を見つけました. `host_to_guest` のエントリが `-1` の場合,対応するホスト頂点にはゲスト頂点が写像されていないことを意味します. ## `einsum` を使った簡潔な目的関数 `objective` を構築する二重の辺ループは,各グラフを辺リストではなく **二値の隣接行列**として表現することで,[`qbpp.einsum`](../../advanced/Einstein-sum.md) を使って 書き換えられます.次のように上三角の隣接行列を用意します: - `A_G[a, b] = 1` ($\{a, b\} \in E_G$ のとき,無向ゲスト辺ごとに 1 エントリ) - `A_H[c, d] = 1` ($\{c, d\} \in E_H$ のとき,無向ホスト辺ごとに 1 エントリ) すると目的関数 $\sum_{(u_G,v_G)\in E_G}\sum_{(u_H,v_H)\in E_H} (x_{u_G,u_H}x_{v_G,v_H}+x_{u_G,v_H}x_{v_G,u_H})$ は 2 つの `einsum` 呼び出しの和として書けます: ```{include} /../programFiles/markDown/example/graph/subgraph-isomorphism.md :start-after: :end-before: ``` subscript `"ab,cd,ac,bd->"` はそのまま $\sum_{a,b,c,d} A_G[a,b]\, A_H[c,d]\, x_{a,c}\, x_{b,d}$ を表しており, [`einsum` ドキュメント](../../advanced/Einstein-sum.md)の QAP 形のテンソル縮約と同じパターンです. 2 つ目の呼び出しはホスト軸を入れ替えた $(u_G, v_G) \mapsto (v_H, u_H)$ の対称写像をカバーします （`ac,bd` の代わりに `ad,bc`）. 得られる QUBO 式の整理後の項集合は for ループ版と完全に同じですが, 構築コードは大幅に短くなり,処理が C++ バックエンド内でマルチスレッド実行 されるため,for ループ版で 1 反復ごとに発生する Python の `ctypes` オーバーヘッドを回避でき大幅に高速になります.トレードオフはメモリ使用量 で,辺リスト方式は $|E_G|+|E_H|$ に比例するのに対し,隣接行列表現は $\Theta(M^2 + N^2)$ となるため,非常に疎で巨大なグラフでは for ループ版 の方が依然として有利です. ## matplotlibによる可視化 以下のコードは,ホストグラフ上で部分グラフ同型の解を可視化します: ```{literalinclude} /../programFiles/pythonPrograms/example/graph/subgraph-isomorphism-program2.py :language: python :caption: subgraph-isomorphism-program2.py ``` 写像されたホスト頂点は赤色で表示され,ゲスト辺に対応する辺が強調表示されます. :::