Help us understand the problem. What is going on with this article?

三角形分割の数え上げとランダムサンプリング

本記事は「データ構造とアルゴリズム Advent Calendar 2019」2 日目の記事です.
1 日目は @deltam さんによる「2つの操作のみで全順列を列挙する:対称群のグラフ上のハミルトン路にもとづく順列生成の紹介と実装」です.
3 日目は @921603 さんによる「Proximity search:列挙アルゴリズムの新たな構築手法」です.

列挙と列挙に挟まれた本記事は,列挙と関連の深い数え上げとランダムサンプリングを扱います.
また,本記事で扱う三角形分割は, 19 日目の @tsukasa__diary さんによる「ドロネー三角形分割の期待最速アルゴリズム(仮)」にも登場すると思います.

はじめに

平面上の $n$ 点の集合 $S$ が与えられたとき,$S$ の三角形分割とは,$S$ の点どうしをつなぐ線分の集合であって,
どの 2 つの線分も互いに交差しないようなもののうち極大なものです.
このような線分の集合は,平面をいくつかの三角形と 1 つの無限に広い外面に分割します.
おそらくこのことから,三角形分割という名前がついています.

本記事では,以下の問題を考えます.

入力として,平面上の $n$ 点の集合 $S$ が与えられる.
ただし,$S$ のどの 3 点も同一直線上にないとする.このとき,
- 数え上げ:$S$ 上の三角形分割の個数 $\mathsf{tr}(S)$ を求めよ.
- ランダムサンプリング:$S$ 上の三角形分割のうち1つを一様ランダムに出力せよ.

ランダムサンプリングは,テストケースや反例の生成において重要です.
また,数え上げはランダムサンプリングと密接な関係をもちます.

三角形分割の数え上げやランダムサンプリングを行うための愚直な方法は,
すべての三角形分割を列挙することです.
すべて列挙してしまえば,個数はわかるし,ランダムに 1 つ選ぶこともできます.
しかし三角形分割は(どの 3 点も一直線上にないような)どんな $S$ に対しても $\Omega(2.43^n)$ 個存在することが知られている1 ので,すべての三角形分割を列挙すると必ず $\Omega(2.43^n)$ 時間かかります.
では,すべての三角形分割を列挙するよりも高速に,数え上げやランダムサンプリングができるでしょうか?

答えは Yes です.

$\textsf{tr}(S)$ を求める時間・空間計算量 $\mathrm{O}(2^n n^2)$ のアルゴリズムが存在する.2
また,同じ時間・空間計算量の前処理ののち,三角形分割を 1 つあたり $\mathrm{O}(n^2)$ 時間でランダムサンプリングできる.

先に述べたとおり,任意の $n$ 点の集合に対して,三角形分割の個数は $\Omega(2.43^n)$ であることが知られています.
したがって,上記のアルゴリズムは,三角形分割すべてを列挙するどんなアルゴリズムよりも指数的に高速であると言えます.

このアルゴリズムは Alvarez and Seidel によって 2013 年の Symposium on Computational Geometry (SoCG) という
国際会議で発表され,その年の best paper award を受賞しました.3
すべての三角形分割を列挙するよりも真に高速に数え上げ・ランダムサンプリングができる
初めてのアルゴリズムであるという点が評価されたようです.
本記事では彼らのアルゴリズムを紹介したいと思います.

注:数え上げに関しては,より高速な $2^{\mathrm{O}(\sqrt{n}\log n)} = 2^{o(n)}$ 時間のアルゴリズムが存在します.4
ただし,そちらはランダムサンプリングに使えるかは(少なくとも本記事の筆者には)明らかではないのと,
full version は約 50 ページあって読むのが大変なので,本記事では Alvarez and Seidel の結果を紹介することにします.

アルゴリズム

$S = \{p_1, \dots, p_n\}$ を平面上の $n$ 点の集合とします.
ただし $S$ に属するどの 3 点も同一直線上にないとします.
さらに,どの 2 点も $x$ 座標を共有しないとします.
このとき,一般性を失わず,$S$ の要素は $x$ 座標の昇順に番号付けられていると仮定できます.
つまり,添字が小さい点ほど左にあるとします.
$S$ の monotone chain とは,パスをなす非交差な線分の集合であって,
最左の点 $p_1$ と最右の点 $p_n$ とを端点とし,どの垂直線にも高々 1 回交わるようなものです.
(直感的には,$p_1$ から $p_n$ に向けて右に進むだけで,左に一度も戻らないパスです.)

$T$ を $S$ の三角形分割とします.
$T$ の monotone chain とは,$S$ の monotone chain であって,その線分がすべて $T$ の辺であるようなものです.

$T$ の monotone chain $C$ に対し, $\Delta_T(C)$ を $C$ の下側にある三角形の集合とします.
ここで「下側」とは,三角形の内部の点すべてが $C$ の下側にあることをいいます.
$T$ の三角形 $t$ が $C$ の advance であるとは, $t$ が $C$ の上側にあり,かつ
$T$ のある monotone chain $C'$ について $\Delta_{T}(C) \cup \{t\} = \Delta_{T}(C')$ がなりたつことを言います.
このとき,$C$ は $C'$ に進むといいます.
三角形 $t$ の頂点を $p_i, p_j, p_k$ とします.ただし,$i < j < k$($p_i, p_j, p_k$ はこの順で左から並んでいる)とします.
このとき, $p_j$ は線分 $[p_i, p_k]$ の上側にある場合と下側にある場合の 2 通りにわかれます.
それぞれの場合について,本記事では $t$ を Type 1 advance, Type 2 advance と呼ぶことにします.
Type 1 のとき,$t$ は $C$ と線分 $[p_i, p_k]$ を共有し,
$C'$ は $C$ から線分 $[p_i, p_k]$ を削除し新たに 2 つの線分 $[p_i, p_j], [p_j, p_k]$ を追加することで得られます.
Type 2 のとき,$t$ は $C$ と線分 $[p_i, p_j], [p_j, p_k]$ を共有し,$C'$ は $C$ からそれらの線分を削除し
新たに線分 $[p_i, p_k]$ を追加することで得られます.

例えば,次の図の三角形分割 $T$ において,赤線で示した線分の集合は $T$ の monotone chain $C$ をなします.
$C$ に対し,青の三角形は Type 1 advance, 緑の三角形は Type 2 advance です.

advance.png
▲ 図 1:monotone chain と advance の例

advance について,以下の補題が成り立ちます.

補題 1:$T$ を $S$ のある三角形分割とする.$T$ の任意の monotone chain $C$ に対し,
$T$ が $S$ の凸包の上側部分でない限り,advance が存在する.

ここで凸包の上側(または下側)部分とは,凸包を最左の点 $p_1$ と最右の点 $p_n$ でわけてできる 2 つの線分の列のうち上側(または下側)にあるもののことを言います.

補題 1 より,$T$ の任意の monotone chain $C$(凸包の上側部分となるものを除く)に対し,advance が存在します.
一般に advance は複数存在しますが,そのうち最も左にあるものは一意に定まります.
したがって,任意の三角形分割 $T$ に対し,以下を満たす $T$ の monotone chain の列 $C_0, \dots, C_M$ が一意に存在します:
1. $C_0$ は凸包の下側部分
2. $C_M$ は凸包の上側部分
3. 各 $C_{\ell}$ は $C_{\ell - 1}$ を一番左の advance によって進めたもの
ここで,$M$ は $T$ の三角形の数です.
($n$ を $S$ の点の個数,$h$ を凸包に含まれる点の個数とすると,$M = 2n - h - 2$ が成り立ちます.)
上記の monotone chain の列を $T$ の leftmost advancing sweep と呼ぶことにします.

次の図に,三角形分割に対する leftmost advancing sweep の例を示します.
先の図と同様に,赤い線分の集合が monotone chain を表します.
図の一番下が $C_0$, 一番上が $C_M = C_4$ です.
また,青の三角形は Type 1 advance を,緑の三角形は Type 2 advance を表します.

leftmost_advancing_sweep.png
▲図 2:leftmost advancing sweep の例

leftmost advancing sweep は $S$ の三角形分割と1対1に対応します.
したがって,$\mathsf{tr}(S)$ を求めるには,leftmost advancing sweep の数を数えれば十分です.
そのために我々は,有向非巡回グラフ(DAG)5 $G_S$ を構築します.
ここで,$G_S$ は source(入次数 0 のノード)と sink(出次数 0 のノード)を 1 つずつもち,
source から sink へのパスと leftmost advancing sweep とが 1 対 1 に対応します.
$G_S$ を構築すると,その source から sink へのパスを数え上げることで leftmost advancing sweep の数,
したがって $\textsf{tr}(S)$ を求めることができます.
また,$G_S$ を用いて効率よくランダムサンプリングを行うことも可能です.

DAG $G_S$ の具体的な定義に移りましょう.
まず,頂点を定義します.
$G_S$ の各頂点は,$S$ に対する marked monotone chain と対応します.
$S$ に対する marked monotone chain は対 $(C, \ell)$ で表されます.
ここで $C$ は $S$ に対する monotone chain です.
また, $\ell$ は整数で,$C$ の左から $\ell$ 番目の線分 $e_\ell$ に印がついていることを表します.
mark $\ell$ は,"leftmost" advancing sweep であることを担保するために管理します.

次に $G_S$ の辺を定義します.
$G_S$ の辺は,ある marked monotone chain が,ある advance によって,
次の marked monotone chain に進むことを表します.
後者の marked monotone chain を,前者の successor といいます.
advance には Type 1 と Type 2 の 2 つがありました.
まずは Type 1 advance に対応する successor の関係を定めます.
$(C, \ell)$ を marked monotone chain とします.
整数 $m \ge \ell$ に対し,$e_m = [p_i, p_j]$ を $C$ の左から $m$ 番目の線分とします.
$S$ の点 $p_{\mu}$ が線分 $[p_i, p_j]$ の"間"にあり($i < \mu < j$),かつ上側にあるとします.
さらに,$p_i, p_{\mu}, p_j$ を頂点とする三角形は内部に $S$ の点を含まないとします.
このとき,$C$ から $e_m$ を削除し,線分 $[p_i, p_{\mu}], [p_{\mu}, p_j]$ を追加して得られる monotone chain を $C'$ とします.
このようにして得られる $(C', m)$ を $(C, \ell)$ の successor とします.

次に, Type 2 advance に対応する successor を定めます.
整数 $m \ge \ell$ に対し,$C$ が 2 つの隣り合う線分 $e_{m-1}, e_m$ を持ち,
$e_{m-1}, e_m$ が張る三角形が $C$ の上側にあり,かつその三角形が内部に $S$ の点を含まないとします.
このとき,$C$ から線分 $e_{m-1}, e_m$ を削除し,$e_{m-1}$ の左の端点と $e_m$ の右の端点とをつなぐ線分を追加して得られる monotone chain を $C'$ とします.
$(C, \ell)$ に対し,$(C', m)$ を successor とします.

DAG $G_S$ は,$(C', m)$ が $(C, \ell)$ の successor であるとき,またそのときのみ $(C, \ell)$ から $(C', m)$ への辺をもちます.

$B$ を $S$ の一番下の monotone chain(つまり,凸包の下側部分)とします.
同様に,$U$ を $S$ の一番上の monotone chain(つまり,凸包の上側部分)とします.
$G_S$ の source を $(B, 1)$ とすると, $(B, 1)$ から $(U, k)$ ($k$ は整数)へのパスが,
$S$ のある三角形分割の leftmost advancing sweep と 1 対 1 に対応します.
便宜上,$G_S$ に頂点 $\top$ を追加し,各 $(U, k)$ から $\top$ に辺を追加します.
また,$\bot := (B, 1)$ とします.
すると,$G_S$ において,$\bot$-$\top$ パスが $S$ のある三角形分割の leftmost advancing sweep と 1 対 1 に対応します.

DAG の例を次の図に示します.
この図では,黒線と赤線で描かれた線分の集合が marked monotone chain を表します.
色を無視した線分の集合が monotone chain であり,赤い線分が mark を表します.
DAGにおいて,$\bot$-$\top$ パスが,leftmost advancing sweep と 1 対 1 に対応していることが確かめられると思います.
例えば一番左の $\bot$-$\top$ パスは図 2 の三角形分割に対する leftmost advancing sweep と対応しています.

DAG.png
▲図 3:DAG $G_S$ の例

詳しいことは後述しますが,この DAG を使って,効率よく数え上げとランダムサンプリングができます.
三角形分割の個数を数え上げるには,DAG 上の $\bot$-$\top$ パスを数え上げればよいです.
これは簡単な動的計画法を使って DAG のサイズに対して線形時間でできます.
また,数え上げの情報を使うとランダムサンプリングも効率よく行えます.

では,$G_S$ のサイズと構築するための計算量はどうなるでしょうか?
次の定理が成り立ちます.

定理 2:DAG $G_S$ のサイズは $\mathrm{O}(2^n n^2)$ である.
また,同じ時間・空間計算量で $G_S$ を構築することができる.

証明:
$G_S$ のサイズを見積もるために,marked monotone chain の種類数を解析します.
$S$ の monotone chain は chain が通る点の集合,すなわち $S$ の部分集合と 1 対 1 に対応するため,その数は高々 $2^n$ です.
また,mark は高々 $n-1$個です.したがって,$G_S$ の頂点数は $\mathrm{O}(2^n n)$ です.
また,各頂点(marked monotone chain)は高々 $n$ 個の successor をもつので,$G_S$ の辺数は $\mathrm{O}(2^n n^2)$ です.
したがって,$G_S$ のサイズは $\mathrm{O}(2^n n^2)$ です.
計算量については,各三角形が内部に点を含むかを多項式時間で前計算しておくと,ある marked monotone chain の successor すべてを $\mathrm{O}(n)$ 時間で生成できます.
したがって,$G_S$ を $\mathrm{O}(2^n n^2)$ 時間・空間で構築することができます.(証明終)

DAG を用いた数え上げとランダムサンプリング

DAG $G_S$ を用いて,数え上げとランダムサンプリングができます.

まず数え上げを考えます.
$G_S$ の頂点 $\alpha$ に対し,$c(\alpha)$ を $\alpha$ から $\top$ へのパスの個数とします.
ただし, $c(\top) = 1$ とします.
$\alpha \neq \top$ のときは $c(\alpha) = \sum_{(\alpha, \beta) \in E} c(\beta)$($E$ は $G_S$ の辺集合)と計算できます.
この式を使うと,各頂点 $\alpha$ を逆トポロジカル順に訪れることで,すべての $c(\alpha)$ を計算できます.
三角形分割の個数は $c(\bot)$ に一致するので,これで三角形分割の個数を求めることができます.
計算時間は $G_S$ のサイズに対して線形時間です.

数え上げの情報を使うと,ランダムサンプリングも効率よくできます.
まず,1 以上 $\mathsf{tr}(S) = c(\bot)$ 以下の整数 $r$ をランダムに生成します.
次に,$G_S$ の「$r$ 番目の」$\bot$-$\top$ パスを見つけます.
このためには, $\bot$ から始めて,以下のように $G_S$ を走査すればよいです.
頂点 $\alpha$ にいるとき,$\alpha$ の successor たちに適当に順序を付け,$\beta_1, \dots, \beta_k$ とします.
$i = 1$ で初期化し,$r > c(\beta_i)$ である限り $r \gets r - c(\beta_i),\, i \gets i + 1$ とします.
$r \le c(\beta_i)$ となった最初の $\beta_i$ に進みます.
この手続きを繰り返すことで,$G_S$ の 「$r$ 番目の」$\bot$-$\top$ パスを見つけることができます.
$r$ は 1 以上 $\mathsf{tr}(S)$ 以下のランダムに生成された整数なので,
これで三角形分割を一様ランダムに生成することができます.
計算時間は,$\bot$-$\top$ パスに含まれる頂点数が $\mathrm{O}(n)$, 各頂点の出次数が $\mathrm{O}(n)$ なので,$\mathrm{O}(n^2)$ です.

おわりに

本記事では,与えられた点集合に対し,三角形分割の数え上げ・ランダムサンプリングを行うアルゴリズムを紹介しました.
アルゴリズムは比較的シンプルなので実装してみたところ,$n = 20$ 個のランダムな点配置に対し,3 秒程度で三角形分割の個数を求めることができました.(2575703778 個,2699808226 個,1569415331 個...など)
使用した計算機は MacBook Pro, 2.3 GHz Intel Core i5, RAM 8 GB です.
計算量 $\mathrm{O}(2^n n^2)$ を考えると,まあ妥当な性能だと思います.
論文では効率良い枝刈りのテクニックについても書かれているのですが,それはまだ本記事の筆者は実装していないので,
それを実装するともっと速くなるかもしれません.
論文の実験結果では $n = 50$ 個のランダムな点配置でも, 19 時間程度かければ数え上げができたと報告されています.

補足として,本記事で紹介したアルゴリズムは三角形分割に特化したものですが, 2014 年に Wettstein によって,三角形分割以外の幾何オブジェクトにも適用できる形に一般化されています.6
また,はじめにの注でも書きましたが,2016 年には Marx and Miltzow によって,より高速な $2^{\mathrm{O}(\sqrt{n}\log n)} = 2^{o(n)}$ 時間のアルゴリズムが提案されています.4

最後までお読みいただき,ありがとうございました.


  1. Micha Sharir, Adam Sheffer, and Emo Welzl. "On degrees in random triangulations of point sets." Journal of Combinatorial Theory, Series A, 118(7):1979-1999, 2011. https://doi.org/10.1016/j.jcta.2011.04.002

  2. 数え上げだけなら空間計算量は $\mathrm{O}(2^n n)$ に削減できますが,本記事では説明を割愛します. 

  3. Victor Alvarez and Raimund Seidel. "A simple aggregative algorithm for counting triangulations of planar point sets and related problems." SoCG 2013. http://doi.org/10.1145/2462356.2462392 

  4. Dániel Marx and Tillmann Miltzow. "Peeling and nibbling the cactus: Subexponential-time algorithms for counting triangulations and related problems." SoCG 2016. http://doi.org/10.4230/LIPIcs.SoCG.2016.52 

  5. 有向非巡回グラフ.つまり,有向閉路を含まない有向グラフ. 

  6. Manuel Wettstein. "Counting and enumerating crossing-free geometric graphs." SoCG 2014. http://doi.org/10.1145/2582112.2582145 

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away