最大マッチング問題¶
無向グラフにおけるマッチングとは,どの2つの辺も共通のノードを持たないような辺の集合です. 無向グラフ \(G=(V,E)\) が与えられたとき,最大マッチング問題は,辺の数が最大となるマッチング \(S \subseteq E\) を求めることを目的とします.
グラフが \(n\) 個の頂点と \(m\) 本の辺を持ち,辺に \(0,1,\ldots,m-1\) のラベルが付いているとします. \(m\) 個のバイナリ変数 \(x_0, x_1, \ldots, x_{m-1}\) を導入し,\(x_i=1\) は辺 \(i\) が選択されている（つまり \(S\) に属する）ことを表します (\(0\le i\le m-1\)). 目的は,選択された辺の数を最大化することです：
マッチング条件を強制するために,共通のノードを持つ選択された辺のペアにペナルティを課します. \(\mathcal{P}\) を,共通の端点を持つ異なる辺の順序なしペア \((e_1,e_2)\) の集合とします. すると,以下のペナルティは,選択された辺がマッチングを形成する場合にのみ値 \(0\) をとります：
目的関数とペナルティを組み合わせて QUBO 式 \(f\) を構築します：
ここで,ペナルティ項に 2 を掛けることで,マッチング制約の違反が目的関数の増加よりもコストが高くなるようにしています. \(f\) を最小化する割り当ては,\(G\) の最大マッチングに対応します.
最大マッチングの Hi-QUBO プログラム
上記の定式化に基づき,以下の Hi-QUBO プログラムは 16 ノードのグラフに対する QUBO 式 \(f\) を構築し,Exhaustive Solver を用いて解きます.
#include <qbpp/qbpp.hpp>
#include <qbpp/exhaustive_solver.hpp>
#include <qbpp/graph.hpp>
int main() {
const size_t N = 16;
std::vector<std::pair<size_t, size_t>> edges = {
{0, 1}, {0, 2}, {1, 3}, {1, 4}, {2, 5}, {2, 6}, {3, 7},
{3, 13}, {4, 6}, {4, 7}, {4, 14}, {5, 8}, {6, 8}, {6, 12},
{6, 14}, {7, 14}, {8, 9}, {9, 10}, {9, 12}, {10, 11}, {10, 12},
{11, 13}, {11, 15}, {12, 14}, {12, 15}, {13, 15}, {14, 15}};
const size_t M = edges.size();
auto x = qbpp::var("x", M);
auto objective = qbpp::sum(x);
auto constraint = qbpp::toExpr(0);
for (size_t i = 0; i < M; ++i) {
for (size_t j = i + 1; j < M; ++j) {
if (edges[i].first == edges[j].first ||
edges[i].first == edges[j].second ||
edges[i].second == edges[j].first ||
edges[i].second == edges[j].second) {
constraint += x[i] * x[j];
}
}
}
auto f = -objective + 2 * constraint;
f.simplify_as_binary();
auto solver = qbpp::ExhaustiveSolver(f);
auto sol = solver.search();
std::cout << "objective = " << objective(sol) << std::endl;
qbpp::graph::GraphDrawer graph;
for (size_t i = 0; i < N; ++i) {
graph.add_node(qbpp::graph::Node(i));
}
for (size_t i = 0; i < M; ++i) {
auto edge = qbpp::graph::Edge(edges[i].first, edges[i].second);
if (sol(x[i])) {
edge.color(1).penwidth(2.0);
}
graph.add_edge(edge);
}
graph.write("maxmatching.svg");
}
このプログラムは式 objective,constraint,f を作成します.f は objective の符号反転にペナルティ項を加えたものです.
Exhaustive Solver が f を最小化し,最適な割り当てが sol に格納されます.
解を可視化するために,GraphDrawer オブジェクト graph を作成し,ノードと辺を追加します.
この可視化では,\(S\) に属する選択された辺（つまり \(x_i=1\) の辺 \(i\)）が強調表示されます.
結果のグラフは描画され,ファイル maxmatching.svg に保存されます：
無向グラフにおけるマッチングとは,共通のノードを持たない辺の集合のことです. 無向グラフ \(G=(V,E)\) が与えられたとき,最大マッチング問題は,辺の数が最大となるマッチング \(S \subseteq E\) を求めることを目的とします.
グラフが \(n\) 個の頂点と \(m\) 本の辺を持ち,辺は \(0,1,\ldots,m-1\) とラベル付けされているとします. \(m\) 個のバイナリ変数 \(x_0, x_1, \ldots, x_{m-1}\) を導入し,\(x_i=1\) は辺 \(i\) が選択されている（すなわち \(S\) に属する）場合にのみ成り立ちます（\(0\le i\le m-1\)）. 目的は選択された辺の数を最大化することです：
マッチング条件を課すために,共通のノードを持つ選択された辺の対にペナルティを与えます. \(\mathcal{P}\) を,共通の端点を持つ異なる辺の順序なし対 \((e_1,e_2)\) の集合とします. すると,以下のペナルティは,選択された辺がマッチングを形成する場合にのみ \(0\) をとります：
目的関数とペナルティを以下のように組み合わせて QUBO 式 \(f\) を構築します：
最大マッチングの PyQBPP プログラム
import pyqbpp as qbpp
N = 16
edges = [
(0, 1), (0, 2), (1, 3), (1, 4), (2, 5), (2, 6), (3, 7),
(3, 13), (4, 6), (4, 7), (4, 14), (5, 8), (6, 8), (6, 12),
(6, 14), (7, 14), (8, 9), (9, 10), (9, 12), (10, 11),(10, 12),
(11, 13),(11, 15),(12, 14),(12, 15),(13, 15),(14, 15)]
M = len(edges)
x = qbpp.var("x", shape=M)
objective = qbpp.sum(x)
constraint = 0
for i in range(M):
for j in range(i + 1, M):
if (edges[i][0] == edges[j][0] or edges[i][0] == edges[j][1] or
edges[i][1] == edges[j][0] or edges[i][1] == edges[j][1]):
constraint += x[i] * x[j]
f = -objective + 2 * constraint
f.simplify_as_binary()
solver = qbpp.ExhaustiveSolver(f)
sol = solver.search()
print(f"objective = {sol(objective)}")
print("Selected edges:", end="")
for i in range(M):
if sol(x[i]) == 1:
print(f" {edges[i]}", end="")
print()
このプログラムは式 objective,constraint,f を作成します.f は objective の符号反転にペナルティ項を加えたものです.
Exhaustive Solver が f を最小化し,最適な割り当てが sol に格納されます.
matplotlib による可視化
以下のコードは,最大マッチングの解を可視化します：
import matplotlib.pyplot as plt
import networkx as nx
G = nx.Graph()
G.add_nodes_from(range(N))
G.add_edges_from(edges)
pos = nx.spring_layout(G, seed=42)
edge_colors = ["#e74c3c" if sol(x[i]) == 1 else "#cccccc"
for i in range(M)]
edge_widths = [2.5 if sol(x[i]) == 1 else 1.0 for i in range(M)]
nx.draw(G, pos, with_labels=True, node_color="#d5dbdb", node_size=400,
font_size=9, edge_color=edge_colors, width=edge_widths)
plt.title("Maximum Matching")
plt.savefig("max_matching.png", dpi=150, bbox_inches="tight")
plt.show()
選択されたマッチング辺は赤で表示されます.