Graphillion で半順序の個数を数えてみた話

Graphillion で遊んでみた。

やったこと

有限集合上の半順序の個数を数えてみた(cf. OEIS A001035)。これを効率的に求めるのは未解決問題らしい。ちなみに有限集合上の位相の個数を数える問題とだいたい同じ。

素数 n を固定して,Graphillion を使って [1..n] 上の半順序をわりと愚直に列挙した。

実行時間は n=6 で 0.3s, n=7 で 8.6s くらい。n=8 になると三時間ぐらい走らせても終わらなかったのでやめた。ちなみに OEIS を見る限り少なくとも n=18 までは既知らしいので,こういうアプローチで新しい値を見つけるのは無理がありそう。

やりかた

a <= b だったら a から b へ辺を伸ばすことにすると,半順序は有向グラフとみなせる。そこで n 点からなる完全グラフの部分グラフで,半順序になっているようなものを列挙すればよい。Graphillion は基本的な集合演算をサポートしているので,基本的には順序の条件をそのまま書くだけでよい。

ただし,Graphillion は無向グラフを扱うライブラリなのに対して,ここでは有向グラフを扱いたい。そこで,source と target を別の点にして有向グラフをむりやり無向グラフで表すことにした。

一応ちゃんと書くと以下のとおり。有向グラフ  (G, E) に対して無向グラフ  (G+G', E') (a, b') \in E' \iff (a, b) \in E となるように作る。ただし  G' G のコピーで, (-)'\colon G \to G'全単射とする。

実際のコードでは,G の点を正の整数,G' の点を負の整数で表した。あと半順序は反射的反対称的推移的関係だけど,対角線部分はあっても情報が増えないので,非対称的推移的関係(strict order; 日本語だと強順序?)にした。こっちのほうが辺の数が減って少し効率がよくなるようだ。

from graphillion import GraphSet

def posets(n):
    univ = [(i, j) for i in range(1, n+1) for j in range(-n, 0)]
    GraphSet.set_universe(univ)
    x = GraphSet({})

    # asymmetry
    for i in range(1, n+1):
        x = x.excluding((i, -i))

    # transitivity
    for i in range(1, n+1):
        for j in range(1, n+1):
            if i != j:
                for k in range(1, n+1):
                    if j != k:
                        y = x.including([(i, -j), (j, -k)]).excluding((i, -k))
                        x = x - y

    return x

import sys
print(len(posets(int(sys.argv[1]))))

実行してみた。

$ time python3 posets.py 6
130023

real	0m0.328s
user	0m0.293s
sys	0m0.028s

$ time python3 posets.py 7
6129859

real	0m8.648s
user	0m8.251s
sys	0m0.324s

出てきた数字は合ってる。

感想とか

ZDD を使うとコンパクトに表現できるかと思ったが,期待したほどコンパクトにはならなかったようだ。ZDD で効率的に表現できるのは基本的には疎な組合せ集合(正確にはちょっと違うと思う)なのだけど,考えてみれば順序全体ってそれほど疎ではないかもしれない。

強引に無向グラフにしてるところにちょっと無駄がある感じはして,Graphillion を使わないで ZDD を自分で操作すればもうちょっと効率よく列挙できそうには思う。

あとは,順序の表現を単に strict order にしたけど,推移律から推論できる辺を除くなどすればもっと辺の数を節約することはできるかもしれない。あるいは,実は BDD のほうがいいとか,(Z)SDD を使うといいとかはあるかもしれない。