概要
「頂点数がNで形が異なるグラフはいくつあるか？」という自然で素朴な疑問について考えます。
本記事は4つのパートに分かれています。
まず、この疑問をラベルなしグラフの数え上げとして定義します。
次に、定義通りに全探索するプログラムとその改善を検討します。
そのあと、上界と下界について考察し、その精度を実験的に確かめます。
最後に、厳密解を（頂点数$N$の分割数を$p(N)$として）$\mathcal{O}(p(N)N)$で求めるアルゴリズムの存在を紹介をします。ただしそのための高度な考察と証明は本記事では行いません。
前提
グラフ理論、集合論に関する初歩的な用語や記法を説明なく用いることがあります。群論の知識はあまり仮定しません。
本記事で「グラフ」とは単純無向グラフを指します。非連結グラフを含みます。
本記事で「頂点数$N$のグラフ」とは、頂点集合が$\lbrace1, 2, ..., N\rbrace$のグラフを指します（つまり、頂点数$N$の任意のグラフの頂点集合が等しいです）
同一なグラフ
2つのグラフ$G, G'$が全く同じ頂点集合・辺集合を持っている時、その2つのグラフは同一であると考えます。辺の曲がり方が違う、頂点の距離が違うなどの違いは無視します。
例えば以下の2つのグラフは同一です。
ラベル付きグラフ (labeled graph)
頂点にラベルが与えられているグラフを指します。
2つのグラフが同一でない場合に区別して数えます。
例えば頂点数3のラベル付きグラフは上の8個です。
同型なグラフ
2つのグラフ$G=\lbrace V, E\rbrace$と$G'=\lbrace V', E'\rbrace$について、以下の条件を満たす全単射写像$f:V \to V'$が存在する時、その2つのグラフは同型であるといいます。
\forall x, y \in V \ (\lbrace x, y\rbrace \in E \Leftrightarrow \lbrace f(x), f(y)\rbrace \in E')
「$G'$は$G$の同型グラフである」というような使い方もします。
また、本記事で問題にするときは常に$V=V'$です。
例えば以下の2つのグラフは同一ではないですが同型です。
自然な定義で、よく分からない人も頂点数3のラベル付きグラフとラベルなしグラフの例を見比べれば納得できると思います。
ラベルなしグラフ
頂点にラベルが与えられていないグラフを指します。
2つのグラフが同型でない場合に区別して数えます。
例えば頂点数3のラベルなしグラフは上の4個です。
頂点数4のラベルなしグラフは上の11個です。
ラベル付きグラフの個数
さて、まずは頂点数が$N$のラベル付きグラフの個数を考えます。
単純無向グラフなので、「異なる2頂点間に辺が0本か1本」あります。異なる2頂点は$\binom{N}{2}$通りあって、それぞれ辺があるかないか2通りなので、頂点数$N$のラベル付きグラフの個数は以下で求まります。
2^\binom{N}{2}
$N=3$で8個、$N=4$で64個、$N=5$で1024個あります。
ラベルなしグラフの個数
本題の、ラベルなしグラフの数え上げについて考えます。
直接閉じた式は思いつきそうにないですね。
定義通りの数え上げ
まずは素直に数え上げるプログラムを書いてみます。
ラベル付きグラフの列挙と、2つのグラフの同型性判定ができればよく、それらを定義通りにプログラムに落とすのは比較的容易です。
ラベル付きグラフの列挙は、$\binom{N}{2}$通りの辺に対してビット全探索すればよいです。
2つのグラフの同型性判定は、$f:V \to V'$を全列挙して定義を満たす写像か判定すれば確実です。ここで考えているのは$V=V'$で、$f:V \to V$の全単射写像とはV上の置換（permutation）なので、順列全探索です。
自明な枝刈りを追加して、以下のようなプログラムで求めることができます。
from itertools import combinations, permutations
from typing import Set, Tuple
def generate_graphs(n):
"""n頂点グラフの辺集合を返すジェネレータ関数
単純無効グラフの個数分生成される（2^(nC2)通り）
"""
possible_edges = list(combinations(range(n), 2))
for i in range(2 ** len(possible_edges)):
# iの2進数表現で1が立っている辺からなるグラフ
edges = set()
for j, edge in enumerate(possible_edges):
# j桁目が立っていたらj番目の辺を追加
if i & (1 << j):
edges.add(edge)
yield edges
def is_isomorphic(
edges1: Set[Tuple[int, int]], edges2: Set[Tuple[int, int]], n: int
) -> bool:
"""グラフが同型か判定する"""
# 辺の本数が異なるグラフは明らかに同型でない
if len(edges1) != len(edges2):
return False
# 同型の定義を満たす全単射写像があるか全探索する
permutations_list = permutations([x for x in range(n)])
for f in permutations_list:
is_currently_isomorphic = True
# 「頂点j, iを結ぶ辺の有無」と「頂点f(j), f(i)を結ぶ辺の有無」が常に一致していればTrue
for i in range(n):
if not is_currently_isomorphic:
break
for j in range(i):
# 実装上j<iなので右辺ではmin,maxを取ればよい
if ((j, i) in edges1) != ((min(f[j], f[i]), max(f[j], f[i])) in edges2):
is_currently_isomorphic = False
break
if is_currently_isomorphic:
return True
return False
def main():
for N in range(8):
start = time.perf_counter()
distinct_graphs = []
for edges1 in generate_graphs(N):
# これまで発見したいずれかのラベルなしグラフと同型な場合はカウントしない
is_new = True
for edges2 in distinct_graphs:
if is_isomorphic(edges1, edges2, N):
is_new = False
break
if is_new:
distinct_graphs.append(edges1)
number_of_distinct_graphs = len(distinct_graphs)
end = time.perf_counter()
print(f"頂点数{N}、ラベルなしグラフの個数：{number_of_distinct_graphs}、実行時間：{end-start:.3f}")
if __name__ == "__main__":
main()
| 頂点数 | ラベルなしグラフの個数 | 実行時間 (秒) |
|---|---|---|
| 0 | 1 | 0.000 |
| 1 | 1 | 0.000 |
| 2 | 2 | 0.000 |
| 3 | 4 | 0.000 |
| 4 | 11 | 0.007 |
| 5 | 34 | 0.057 |
| 6 | 156 | 23.563 |
| 7 | - | （30000秒以上） |
$N=7$は一晩立っても終わりませんでした。
外側のループがラベル付きグラフの個数（$2^\binom{7}{2}=2097152$）回、内側のループがこれまでに発見したラベルなしグラフ回（N=7で最大1044回）、その中で$7!=5040$の置換を調べているので、まぁこれくらいが限界でしょう。
グラフ同型性判定問題（Graph isomorphism problem）
2つのグラフが同型か否かを判定する問題を「グラフ同型性判定問題」と言います。
前節では全ての全単射写像を列挙して判定する実装をしましたが、もう少し賢い方法がありそうですよね。
例えば「辺の本数が異なる場合は明らかに同型でない」としましたが、それ以外にも「次数列が多重集合として異なれば同型でない」というのも明らかです。また、写像をDFSで列挙しながら、途中で明らかに同型でないことが分かった場合（例えば閉路の長さや個数が異なる場合）バックトラックすることで枝刈りするアイディアが考えられます。
あるいは別のアプローチとして、何らかの方法で2つのラベル付きグラフの正規化を定義して、正規化グラフの一致と同型性が同値になるようにして判定するというアイディアも自然です。
調べてみるとどちらもあるようで、前者はVF2++などのアルゴリズム、後者はnautyなどのアルゴリズムが知られています。
いずれも最悪ケースでは指数時間掛かりますが、実用的なグラフに対しては十分高速に動作するとされています。
なお、グラフ同型性判定問題の計算量クラスがPなのかNP完全なのかNP中間（NP-intermediate）なのかというのは未解決問題です。
同型性判定をVF2++に変更
VF2++はNetworkXという、グラフ分野のポピュラーなPythonパッケージに採用されているため、簡単に差し替えられます。
def is_isomorphic(
edges1: Set[Tuple[int, int]], edges2: Set[Tuple[int, int]], n: int
) -> bool:
"""グラフが同型か判定する"""
g1 = nx.Graph()
g1.add_edges_from(edges1)
g2 = nx.Graph()
g2.add_edges_from(edges2)
gm = isomorphism.GraphMatcher(g1, g2)
return gm.is_isomorphic()
すごく簡単ですね（実際は毎回Graphを作るのは無駄なので使いまわせるように、is_isomorphicの引数にGraphを渡すように変更しています）
| 頂点数 | ラベルなしグラフの個数 | 実行時間 (秒) | （参考：愚直実装） |
|---|---|---|---|
| 0 | 1 | 0.000 | 0.000 |
| 1 | 1 | 0.000 | 0.000 |
| 2 | 2 | 0.000 | 0.000 |
| 3 | 4 | 0.001 | 0.000 |
| 4 | 11 | 0.012 | 0.007 |
| 5 | 34 | 0.259 | 0.057 |
| 6 | 156 | 3.233 | 23.563 |
| 7 | 1044 | 1264.953 | （30000秒以上） |
速い！$N=7$が求まるようになりました！！
それでも$N=8$以上は厳しそうですね。$N=8$のラベル付きグラフは$2.7×10^8$通りあり、内側のループも最大$10^4$回回るので、同型性判定が爆速だとしても時間が掛かりすぎます。
「全てのラベル付きグラフを列挙して1つ1つ同型性判定する」という方針の限界でしょう。
ここで厳密解は小休止して、「頂点数$N$のラベルなしグラフは大体どれくらいあるか？」を見積ってみましょう。
ラベルなしグラフの個数の上界
頂点数$N$のラベルなしグラフにとって、頂点数$N$のラベル付きグラフの個数は自明な上界です。
2^\binom{N}{2}
ラベルなしグラフの個数の下界
頂点数$N$のあるグラフ$G$について、$G$と同型なラベル付きグラフが最大いくつ存在するか考えます。
同型の定義より、$G$と同型なラベル付きグラフの個数が最大になるのは、頂点の置換全てが同一でないラベル付きグラフになる場合で、その個数は$N!$個です。
同様に、頂点数$N$の任意のグラフ$G$について同型なラベル付きグラフは$N!$以下なので、頂点数$N$のラベル付きグラフ全体を「同型なグラフ同士を1まとめのグループにする」と、それぞれのグループのサイズは$N!$以下になります。
さて頂点数$N$のラベル付きグラフは$2^\binom{N}{2}$個あり、それぞれのグループのサイズが$N!$以下であることから、グループの個数（＝ラベルなしグラフの個数）は
\frac{2^{\binom{N}{2}}}{N!}
以上であると言えます。
実際の頂点数4のラベルなしグラフと、それと同型なラベル付きグラフの個数です。このそれぞれのグループの個数がいずれも$4!=24個$以下であるということです。
上界、下界、実際の個数の可視化
対数スケールでプロットします。
べらぼうにデカいですね。上界がデカすぎるので麻痺しますが、デカいでお馴染み$2^N$が$N=20$で概ね$10^6$であることを思い出すと、下界もぶっ飛んでます。
ここに実際の個数もプロットしてみましょう。我々はまだ$N\geq8$のラベルなしグラフの厳密な個数を知らないはずですが、ここまでの計算で得た「1,1,2,4,11,34,156,1044」をOEISに打ち込むことで、続きを錬金できます。
import math
import matplotlib.pyplot as plt
# 20頂点以下のラベルなしグラフの個数（https://oeis.org/A000088）
number_of_distinct_graphs = [
1,
1,
2,
4,
11,
34,
156,
1044,
12346,
274668,
12005168,
1018997864,
165091172592,
50502031367952,
29054155657235488,
31426485969804308768,
64001015704527557894928,
245935864153532932683719776,
1787577725145611700547878190848,
24637809253125004524383007491432768,
645490122795799841856164638490742749440,
]
# 上界：2^(nC2)
def calc_upper_bound(n):
return 2 ** math.comb(n, 2)
# 下界：2^(nC2)/n!
def calc_lower_bound(n):
return calc_upper_bound(n) / math.factorial(n)
def main():
MAX_N = 20
assert MAX_N < len(number_of_distinct_graphs)
n = list(range(MAX_N + 1))
ub_list = [calc_upper_bound(x) for x in n]
lb_list = [calc_lower_bound(x) for x in n]
plt.plot(n, ub_list, label="upper bound", color="red", linewidth=1)
plt.plot(n, lb_list, label="lower bound", color="green", linewidth=1)
plt.plot(n, number_of_distinct_graphs[: MAX_N + 1], label="ditsinct graphs")
plt.yscale("log")
plt.title("Number of Distinct Graphs")
plt.legend()
plt.show()
if __name__ == "__main__":
main()
| 頂点数 | 実際の個数 | 下界 | 上界 |
|---|---|---|---|
| $0$ | $1$ | $1.0$ | $1$ |
| $1$ | $1$ | $1.0$ | $1$ |
| $2$ | $2$ | $1.0$ | $2$ |
| $3$ | $4$ | $1.33$ | $8$ |
| $4$ | $11$ | $2.67$ | $64$ |
| $5$ | $34$ | $8.53$ | $1024$ |
| $6$ | $156$ | $45.5$ | $3.28 \times 10^4$ |
| $7$ | $1040$ | $416$ | $2.10 \times 10^6$ |
| $8$ | $1.23 \times 10^4$ | $6658$ | $2.68 \times 10^8$ |
| $9$ | $2.75 \times 10^5$ | $1.89 \times 10^5$ | $6.87 \times 10^{10}$ |
| $10$ | $1.20 \times 10^7$ | $9.70 \times 10^6$ | $3.52 \times 10^{13}$ |
| $11$ | $1.02 \times 10^9$ | $9.03 \times 10^8$ | $3.60 \times 10^{16}$ |
| $12$ | $1.65 \times 10^{11}$ | $1.54 \times 10^{11}$ | $7.38 \times 10^{19}$ |
| $13$ | $5.05 \times 10^{13}$ | $4.85 \times 10^{13}$ | $3.02 \times 10^{23}$ |
| $14$ | $2.91 \times 10^{16}$ | $2.84 \times 10^{16}$ | $2.48 \times 10^{27}$ |
| $15$ | $3.14 \times 10^{19}$ | $3.10 \times 10^{19}$ | $4.06 \times 10^{31}$ |
| $16$ | $6.40 \times 10^{22}$ | $6.35 \times 10^{22}$ | $1.33 \times 10^{36}$ |
| $17$ | $2.46 \times 10^{26}$ | $2.45 \times 10^{26}$ | $8.71 \times 10^{40}$ |
| $18$ | $1.79 \times 10^{30}$ | $1.78 \times 10^{30}$ | $1.14 \times 10^{46}$ |
| $19$ | $2.46 \times 10^{34}$ | $2.46 \times 10^{34}$ | $2.99 \times 10^{52}$ |
| $20$ | $6.45 \times 10^{38}$ | $6.45 \times 10^{38}$ | $1.57 \times 10^{60}$ |
ほぼ実際の個数に一致する良い下界ですね。
Nが大きいほど実際の個数よりも大体の個数が知りたいユースケースが多いと思うので、これでいいんじゃないでしょうか？
めでたしめでたしです。
いや私はどうしても厳密な個数が知りたいです
そんな……どうして……。
ご用意いたしました。
"""
https://oeis.org/A000088
Chai Wah Wu, Jul 02 2024
"""
import time
from itertools import combinations
from math import prod, factorial, gcd
from fractions import Fraction
from sympy.utilities.iterables import partitions
def A000088(n):
return int(
sum(
Fraction(
1
<< sum(p[r] * p[s] * gcd(r, s) for r, s in combinations(p.keys(), 2))
+ sum((q >> 1) * r + (q * r * (r - 1) >> 1) for q, r in p.items()),
prod(q**r * factorial(r) for q, r in p.items()),
)
for p in partitions(n)
)
)
def main():
for i in range(0, 101, 10):
start = time.perf_counter()
print(i, A000088(i))
end = time.perf_counter()
print(f"頂点数{i}の実行時間：{end - start:.3f}")
if __name__ == "__main__":
main()
……？？？？
こちらのコードはOEISから窃盗したものです。
これなら頂点$100$のラベルなしグラフの個数も求まります。
| 頂点数 | ラベルなしグラフの個数 | 実行時間 (秒) |
|---|---|---|
| 0 | 1 | 0.000 |
| 10 | 12005168 | 0.000 |
| 20 | 645490122795799841856164638490742749440 | 0.002 |
| 30 | 334494316309257669249439569928080028956631479935393064329967834887217734534880582749030521599504384 | 0.020 |
| 40 | 7793841167914977954582550817575177766066055272533160501864210580719699592280766598762108507458913936081932965352037372886593259286753883857016383307981863462449691949358853053120648183808 | 0.152 |
| 50 | 18996334831796305871264776333897774457909428769342981426195127107030831238341767361078220789873024149071864408554239370333393376289201491326704031896750415884761063412211252965090574099406295069297305098820585820480974508581202071478261179801957205826115377702867655279140649075352303688667581478179176448 | 0.897 |
| 60 | 7996822858506090743800390365651771608698688069947967132554412354171399322696112483853872208710884748423663796708670753197241378824822980807388078880750196247172406732905727371278913588018333888632848897187108282154641223814874243509088688915549395092057735041789325876615710812148439042842259783997040919940986417289398741439045573791717592119276096665631039828680562427331959523029806766864353627853821416488251739960907314850238929686508536734941184 | 4.663 |
| 70 | 811025469117756850101645165297636968394651070439058650511864910423859549885226389197744046865772849130807350523708001048163948414506716560110394686725323892157505896014142873408147230679071620087896805386165482998746031272578100829012067672123742574426280468994049568047862721186193598987029840188542288417488989225234854755464228426042673199027174195987032080966433597584556769307845316423717670714072507067980701812503312824558796648216394318207710633571527179953231011350496755058729638337084165596553241909397539430082838094716542976099146168079325846625711495090588183972616138993979123159365985809610247764726074218381312 | 21.702 |
| 80 | 25122252462708578937918182926643824995601581876813618155058595516866323395773114190157775032371048608573189104819084848705904868854948334137992363729517482412177300453391604676441864502263055191761799266348963348501714182419601289332061048574704914364838661191748380568338985959818737784559152961655073493809337195872929546048625531837508404652385056319139551869585993645617925054637758412645518666435830113354714296506355716277256535945204807141894964731857335976616713429709547402321240666144124147884061140667056264652766454714026037555679904106808894524715215335318479034572640631212994671350042282534651363833249057171366504990042608461188742439271791862357912725232734924043387792249517070346729165756368170891793343456228230058639080854098686044040705200589614278929650953870987902250164410262886084506783022785194968217550848 | 92.640 |
| 90 | 283920560848950935497722181872546452058706558012643973234739332377178573927628996825196180979821594530457027538738949940778894550094127486078049253021499721100387544954648228890729254583326603955648812896009780171546302360511697528581174356444979068914995016019309074343323782811729043274790596188159542641303746368280313792913285293256635684095220794876652194101554476190705040620928993536578374733810369856886673062337917058403325406678941705449596110350834033559591372890864137449135248041718935175767536514322164691762282767903705536281845521100996995079877810847731312771214343404124140690514191017181638785912870195463668498373645297543926035624939010298061718400648211338841976221130242246565632983095596281147247495445622875367766626744185957010509969459875846313045869586415885340316885615926491134860388229650835388434706047592350863110940206695152288581096084173945794956015624784548453015010284655625529941739680125753532463967359509341608035056518097125411797045433256278613882738527789812813646372996806550004035768496399880303505930867276696577124597760 | 364.663 |
| 100 | 1344234018830050151564780041892128661588217920322175866733047523526572224819608031470152841097730274344158146075356239066674085877666421985093753849051210371487458299332932896874231165042092495362516523377967359346837046787608653334137229336386808737891630136107821433214354213855171893998743007384073524334786906015236575007806511384151113957153639556141939400032090201931586254301901478805218791851843030658508064333241361442393618902028321499459376620207995686971550252753073672772726632339924983691695158972438431612114038475546425749932363148754922032992360984304051927185271025979701973632432290012073991955332121895864434472383448711419186222268140686086115448172198858247173642841694227676231047701650607816758922270034852134640384866512862004699149922492419425425497506922674599831004109661939231353143610118668474253104324772388896830721426097955476885709342337803540179000695713661976647222977100747309218639993723000369469498281782795327711016434608169797477445794481472849129467751866549261655543294615771763917761748371734376410198979492287782333093401525672340085476993472949272255597226688315403931084273461600542204175181025343297470470805066490255888432086427816425649725781478388911692828043651184922785020772722361026974488152599648880251070377167694426556572173316140952186971363878935221832753147690388744568832 | 1336.405 |
def A000088(n):
def calculate_numerator(p):
# 個数が異なるサイクル間の組み合わせ
part1 = 1
for r, s in combinations(p.keys(), 2):
part1 *= 2 ** (p[r] * p[s] * gcd(r, s))
# 個数が同じだけど異なるサイクル間の組み合わせ
part2 = 1
for q, r in p.items():
part2 *= 2 ** (q * comb(r, 2))
# 同じサイクル内の組み合わせ
part3 = 1
for q, r in p.items():
part3 *= 2 ** (q // 2 * r)
# それらの総積
return part1 * part2 * part3
def calculate_denominator(p):
return prod(q**r * factorial(r) for q, r in p.items())
ans = 0
for p in partitions(n):
# 分割のパターンごとに数え上げる
ans += Fraction(calculate_numerator(p), calculate_denominator(p))
return ans
一応コードを読み解くとこんな感じです。
『いかにして問題を解くか』でお馴染みポリア先生の数え上げ定理などを用いた高度な考察が必要で、自分の手には余るため、最後に参考文献案内をするに留めます。
さて、いくつか前計算して高速化できるにしても、頂点数$N$に対して、$N\times (Nの分割数)$くらいは掛かりそうですね。多倍長整数の計算コストも負担になっています。
分割数のOEISによると、$100$の分割数は約$2\times 10^8$、$200$の分割数は約$4\times 10^{12}$なので、あまりにも巨大な頂点に対して厳密解を求めるのは難しそうです。
一方で下界による評価であれば、1000頂点とかでも爆速で求まって嬉しいです。$7.554\times 10^{147796}$個らしいです。嬉しいですね！！
本記事で作成したプログラム
関連、参考文献
ラベルなしグラフの数え上げ
ABC284-Ex - Count Unlabeled Graphs
$K=1$の場合、本記事で扱ったラベルなしグラフの数え上げです。$N$は最大$30$なので最後に述べた考察が必要です。とても難しいです（diff3125）（10000人参加して解けたのは16人）。
解説動画で「異なるサイクル間の組み合わせでgcdが出てくる理由」など視覚的に解説してくれています。
[Educational] Combinatorics Study Notes (4) -Black_Fate's blog(codeforces)
最終的にポリアの数え上げ定理を使ってABC284-Exを解くまでの組み合わせ論勉強ノートです。
『グラフの数え上げ 母関数を礎にして』 共立出版
第4章でラベルなしグラフ（非標識グラフ）の数え上げを扱います。日本語で最も詳しい本だと思います。
いろいろなグラフの数え上げ
例えば「ラベルなし木」ならどうか？など様々なグラフ数え上げがまとまっています。
ラベルなしグラフの個数の下界
『離散数学への招待 上』 丸善出版
下界の議論は本書が元ネタです。「グラフを描画して実際の個数がこの下界に近いことを確かめる」というのも本書の演習問題です。
グラフ同型性判定問題
The Graph Isomorphism Problem -Communications of the ACM
グラフ同型性判定問題の2020年の広いレビュー記事です。
グラフ同型性判定問題への招待 電子情報通信学会
日本語だとこのあたりでしょうか（大学や会社で読める人は）
計算量クラスが未解決問題である点を多めに論じていて、アルゴリズムの概要紹介は少なめです。
VF2++
NetworkXのリファレンス
実装も読めます。
VF2++—An improved subgraph isomorphism algorithm
オリジナルの論文です。
VF2 アルゴリズムの雑メモ
日本語だとこのあたりでしょうか。
Nauty
nauty and Traces
公式です。使い方ドキュメントや提案論文へのリンクも充実しています。