部分グラフ同型問題¶

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)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	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\) 頂点のホストグラフに対する部分グラフ同型問題を解きます:

subgraph-isomorphism-program1.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>
#include <qbpp/graph.hpp>

int main() {
  const size_t N = 10;
  std::vector<std::pair<size_t, size_t>> host = {
      {0, 1}, {0, 2}, {1, 3}, {1, 4}, {1, 6}, {2, 5}, {3, 7}, {4, 6},
      {4, 7}, {5, 6}, {5, 8}, {6, 8}, {6, 7}, {7, 9}, {8, 9}};

  const size_t M = 6;
  std::vector<std::pair<size_t, size_t>> guest = {
      {0, 1}, {0, 2}, {1, 2}, {1, 3}, {2, 3}, {2, 5}, {3, 4}, {4, 5}};

  auto x = qbpp::var("x", M, N);

  auto host_assigned = qbpp::vector_sum(x, 0);

  auto constraint = qbpp::sum(qbpp::vector_sum(x, 1) == 1) +
                    qbpp::sum(0 <= host_assigned <= 1);

  auto objective = qbpp::toExpr(0);

  for (const auto& e_g : guest) {
    for (const auto& e_h : host) {
      objective += x[e_g.first][e_h.first] * x[e_g.second][e_h.second] +
                   x[e_g.first][e_h.second] * x[e_g.second][e_h.first];
    }
  }

  auto f = -objective + constraint * (M * N);

  f.simplify_as_binary();

  auto solver = qbpp::EasySolver(f);
  auto sol = solver.search({{"target_energy", std::to_string(-static_cast<int>(guest.size()))}});

  std::cout << "sol(x) = " << sol(x) << std::endl;

  std::cout << "sol(objective) = " << sol(objective) << std::endl;
  std::cout << "sol(constraint) = " << sol(constraint) << std::endl;

  auto guest_to_host = qbpp::onehot_to_int(sol(x), 1);
  std::cout << "guest_to_host = " << guest_to_host << std::endl;

  auto host_to_guest = qbpp::onehot_to_int(sol(x), 0);
  std::cout << "host_to_guest = " << host_to_guest << std::endl;

  qbpp::graph::GraphDrawer guest_graph;
  for (size_t i = 0; i < M; ++i) {
    guest_graph.add(qbpp::graph::Node(i));
  }
  for (const auto& e : guest) {
    guest_graph.add(qbpp::graph::Edge(e.first, e.second));
  }
  guest_graph.write("guest_graph.svg");

  qbpp::graph::GraphDrawer host_graph;
  for (size_t i = 0; i < N; ++i) {
    host_graph.add(qbpp::graph::Node(i));
  }
  for (const auto& e : host) {
    host_graph.add(qbpp::graph::Edge(e.first, e.second));
  }
  host_graph.write("host_graph.svg");

  qbpp::graph::GraphDrawer graph;
  for (size_t i = 0; i < N; ++i) {
    graph.add(qbpp::graph::Node(i).color(sol(host_assigned[i])));
  }

  std::vector<std::vector<bool>> guest_adj(N, std::vector<bool>(N, false));
  for (auto [u, v] : guest) {
    guest_adj[u][v] = guest_adj[v][u] = true;
  }
  for (const auto& e_h : host) {
    auto u = host_to_guest[e_h.first];
    auto v = host_to_guest[e_h.second];
    if (u != -1 && v != -1 &&
        guest_adj[static_cast<size_t>(u)][static_cast<size_t>(v)]) {
      graph.add(
          qbpp::graph::Edge(e_h.first, e_h.second).color(1).penwidth(2.0f));
    } else {
      graph.add(qbpp::graph::Edge(e_h.first, e_h.second));
    }
  }

  graph.write("subgraph_isomorphism.svg");
}

ゲストグラフとホストグラフは,それぞれ辺リスト 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 に可視化されます.写像で選択されたホスト頂点と,ゲスト辺に対応するホスト辺がハイライトされています.

このプログラムは以下の出力を生成します:

sol(x) = {{0,1,0,0,0,0,0,0,0,0},{0,0,0,0,1,0,0,0,0,0},{0,0,0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0,1,0}}
sol(objective) = 8
sol(constraint) = 0
guest_to_host = {1,4,6,7,9,8}
host_to_guest = {-1,0,-1,-1,1,-1,2,3,5,4}

目的関数値はゲスト辺の数（\(|E_G|=8\)）に等しく,すべての制約が満たされています（constraint = 0）. したがって,プログラムは有効な部分グラフ同型に対応する最適解を見つけました. host_to_guest のエントリが -1 の場合,対応するホスト頂点にはゲスト頂点が写像されていないことを意味します.

einsum を使った簡潔な目的関数

objective を構築する二重の辺ループは,各グラフを辺リストではなく 二値の隣接行列として表現することで,qbpp::einsum を使って書き換えられます.次のように上三角の隣接行列を用意します:

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 呼び出しの和として書けます:

  // 隣接行列を構築（無向辺ごとに 1 エントリ）
  auto A_G = qbpp::array<qbpp::coeff_t>(M, M);
  for (const auto& e : guest) A_G[e.first][e.second] = 1;
  auto A_H = qbpp::array<qbpp::coeff_t>(N, N);
  for (const auto& e : host)  A_H[e.first][e.second] = 1;

  // Σ_{a<b, c<d} A_G[a,b] * A_H[c,d] * (x[a,c]*x[b,d] + x[a,d]*x[b,c])
  auto objective = qbpp::einsum<0>("ab,cd,ac,bd->", A_G, A_H, x, x)
                 + qbpp::einsum<0>("ab,cd,ad,bc->", A_G, A_H, x, x);

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 ドキュメントの 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 ループ版の方が依然として有利です.

例として,以下のホストグラフとゲストグラフを考えます:
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定式化

例えば,部分グラフ同型問題の解は以下の \(6\times 10\) バイナリ行列で表現できます:

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	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} \]

部分グラフ同型問題のPyQBPPプログラム

上記のQUBO定式化に基づき,以下のPyQBPPプログラムは \(M=6\) 頂点のゲストグラフと \(N=10\) 頂点のホストグラフに対する部分グラフ同型問題を解きます:

subgraph-isomorphism-program1.py¶

import pyqbpp as qbpp

N = 10
host = [
    (0, 1), (0, 2), (1, 3), (1, 4), (1, 6), (2, 5), (3, 7), (4, 6),
    (4, 7), (5, 6), (5, 8), (6, 8), (6, 7), (7, 9), (8, 9)]

M = 6
guest = [
    (0, 1), (0, 2), (1, 2), (1, 3), (2, 3), (2, 5), (3, 4), (4, 5)]

x = qbpp.var("x", shape=(M, N))

host_assigned = qbpp.vector_sum(x, axis=0)

constraint = qbpp.sum(qbpp.vector_sum(x, axis=1) == 1) + \
             qbpp.sum((0 <= host_assigned) & (qbpp.same <= 1))

objective = 0
for ug, vg in guest:
    for uh, vh in host:
        objective += x[ug][uh] * x[vg][vh] + x[ug][vh] * x[vg][uh]

f = -objective + constraint * (M * N)
f.simplify_as_binary()

solver = qbpp.EasySolver(f)
sol = solver.search(target_energy=-len(guest))

print(f"sol(x) = {sol(x)}")
print(f"sol(objective) = {sol(objective)}")
print(f"sol(constraint) = {sol(constraint)}")

guest_to_host = qbpp.onehot_to_int(sol(x), axis=1)
print(f"guest_to_host = {guest_to_host}")

host_to_guest = qbpp.onehot_to_int(sol(x), axis=0)
print(f"host_to_guest = {host_to_guest}")

関数 qbpp.onehot_to_int() を用いて,ゲスト頂点からホスト頂点への写像（guest_to_host,\(\sigma\)）とホスト頂点からゲスト頂点への写像（host_to_guest,\(\sigma^{-1}\)）も出力します.

このプログラムは以下の出力を生成します:

sol(x) = {{0,1,0,0,0,0,0,0,0,0},{0,0,0,0,1,0,0,0,0,0},{0,0,0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,1,0,0},{0,0,0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0,1,0}}
sol(objective) = 8
sol(constraint) = 0
guest_to_host = {1,4,6,7,9,8}
host_to_guest = {-1,0,-1,-1,1,-1,2,3,5,4}

目的関数値はゲスト辺の数（\(|E_G|=8\)）に等しく,すべての制約が満たされています（constraint = 0）. したがって,プログラムは有効な部分グラフ同型に対応する最適解を見つけました. host_to_guest のエントリが -1 の場合,対応するホスト頂点にはゲスト頂点が写像されていないことを意味します.

einsum を使った簡潔な目的関数

objective を構築する二重の辺ループは,各グラフを辺リストではなく 二値の隣接行列として表現することで,qbpp.einsum を使って書き換えられます.次のように上三角の隣接行列を用意します:

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 呼び出しの和として書けます:

# 隣接行列を構築（無向辺ごとに 1 エントリ）
A_G = qbpp.array([0] * (M * M), shape=(M, M))
for u, v in guest:
    A_G[u][v] = 1
A_H = qbpp.array([0] * (N * N), shape=(N, N))
for u, v in host:
    A_H[u][v] = 1

# Σ_{a<b, c<d} A_G[a,b] * A_H[c,d] * (x[a,c]*x[b,d] + x[a,d]*x[b,c])
objective = qbpp.einsum("ab,cd,ac,bd->", A_G, A_H, x, x) \
          + qbpp.einsum("ab,cd,ad,bc->", A_G, A_H, x, x)

得られる QUBO 式の整理後の項集合は for ループ版と完全に同じですが, 構築コードは大幅に短くなり,処理が C++ バックエンド内でマルチスレッド実行されるため,for ループ版で 1 反復ごとに発生する Python の ctypes オーバーヘッドを回避でき大幅に高速になります.トレードオフはメモリ使用量で,辺リスト方式は \(|E_G|+|E_H|\) に比例するのに対し,隣接行列表現は \(\Theta(M^2 + N^2)\) となるため,非常に疎で巨大なグラフでは for ループ版の方が依然として有利です.

matplotlibによる可視化

以下のコードは,ホストグラフ上で部分グラフ同型の解を可視化します:

subgraph-isomorphism-program2.py¶

import matplotlib.pyplot as plt
import networkx as nx

G_host = nx.Graph()
G_host.add_nodes_from(range(N))
G_host.add_edges_from(host)
pos = nx.spring_layout(G_host, seed=42)

# 写像されたホスト頂点を特定
mapped = [0] * N
for i in range(M):
    for j in range(N):
        if sol(x[i][j]) == 1:
            mapped[j] = 1
colors = ["#e74c3c" if mapped[j] else "#d5dbdb" for j in range(N)]

# ゲスト辺に対応する辺を強調
guest_adj = set()
for u, v in guest:
    guest_adj.add((u, v))
    guest_adj.add((v, u))

guest_to_host_map = {}
for i in range(M):
    for j in range(N):
        if sol(x[i][j]) == 1:
            guest_to_host_map[i] = j
host_to_guest_map = {v: k for k, v in guest_to_host_map.items()}

edge_colors = []
edge_widths = []
for u, v in host:
    gu = host_to_guest_map.get(u)
    gv = host_to_guest_map.get(v)
    if gu is not None and gv is not None and (gu, gv) in guest_adj:
        edge_colors.append("#e74c3c")
        edge_widths.append(2.5)
    else:
        edge_colors.append("#cccccc")
        edge_widths.append(1.0)

nx.draw(G_host, pos, with_labels=True, node_color=colors, node_size=400,
        font_size=9, edge_color=edge_colors, width=edge_widths)
plt.title("Subgraph Isomorphism")
plt.savefig("subgraph_isomorphism.png", dpi=150, bbox_inches="tight")
plt.show()

写像されたホスト頂点は赤色で表示され,ゲスト辺に対応する辺が強調表示されます.

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0

\(i\)	\(\sigma(i)\)	1	4	6	7	8	9
0	1	1	0	0	0	0	0
1	4	0	1	0	0	0	0
2	6	0	0	1	0	0	0
3	7	0	0	0	1	0	0
4	9	0	0	0	0	0	1
5	8	0	0	0	0	1	0