部分消化問題を解く Go のプログラムを書いてみた話。
あるいは、分岐する再帰関数を Goroutine で並行処理してみた話

なお、私はバイオ系でもアルゴリズム屋でもないので正確性については悪しからず。

部分消化問題とは

整数の集合 $X \subset \mathbb{Z}$ を、多重集合 $\Delta X$ より求めよ。ここで $\Delta X$ は次式で定義される。

\Delta X = \left\{ x - y \mid x, y \in X \land x > y \right\}

例えば、$X = \{0, 2, 4, 7, 10 \}$ なら $\Delta X = \{2, 2, 3, 3, 4, 5, 6, 7, 8, 10\}$である。この、$\{2, 2, 3, 3, 4, 5, 6, 7, 8, 10\}$ から $\{0, 2, 4, 7, 10 \}$ を求めることが、部分消化問題を解くことに当たる。

正格な話を知りたい場合は下記で

バイオインフォマティクスのためのアルゴリズム入門を買う
リンク先を読む

ホモメトリックな解

部分消化問題の解は無数に存在する。そのため、いくつかの解はまとめて一つの解として扱う必要がある。
$\Delta A = \Delta B$ が成り立つとき、 $A$ と $B$ はホモメトリックであると言う。

ずらし

$A$ の全ての要素を一定数 $v$ だけずらした集合 $A \oplus v$ は、 $A$ とホモメトリックである。

$A \oplus \{v\} = \{ a + v \mid a \in A, v \in \{v\}\}$
$\Delta A = \Delta (A \oplus v)$

例えば、$X = \{0, 2, 4, 7, 10 \}$ と $X \oplus \{100\} = \{100, 102, 104, 107, 110\}$ について、$\Delta X = \Delta (X \oplus \{100\}) = \{2, 2, 3, 3, 4, 5, 6, 7, 8, 10\}$ である。

今回の実装では、ずらしによって一致する解はすべて同一視して、最小元がゼロとなるような集合 $A \oplus \{-\min A\}$ のみを解として扱っている。

反射

$A$ の全ての要素の符号を反転させた集合 $- A$ は、 $A$ とホモメトリックである。

$-A = \{ -a \mid a \in A\}$
$\Delta A = \Delta (- A)$

例えば、$X = \{0, 2, 4, 7, 10 \}$ と $- X = \{0, -2, -4, -7, -10\}$ について、$\Delta X = \Delta (-X) = \{2, 2, 3, 3, 4, 5, 6, 7, 8, 10\}$ である。

今回の実装では、基本的に反射は別の解として区別している。すなわち、次の２つを区別している。
ただし、左右対称で、$A \oplus \{-\min A\} = -A \oplus \{\max A\}$となる場合、単一の解としている。

$A \oplus \{-\min A\}$
$-A \oplus \{-\min (-A)\} = -A \oplus \{\max A\}$

例えば、$\{0, 2, 4, 7, 10 \}$ と $\{0, 3, 6, 8, 10 \}$ は区別しているが、$\{0, 2, 5, 8, 10\}$は左右対称であるため、その反射と区別していない（できない）。

自明でない例

任意の２つの集合 $U, V$ について、各集合の要素の和からなる集合 $U \oplus V$ と、各集合の要素の差からなる集合 $U \ominus V$ はホモメトリックである。

$U \oplus V = \{ u + v \mid u \in U, v \in V\}$
$U \ominus V = \{ u - v \mid u \in U, v \in V\}$
$\Delta(U \oplus V) = \Delta(U \ominus V)$

例えば、$U = \{6, 7, 9\}$、$V = \{-6, 2, 6\}$とした場合、次の２つはホモメトリックである。

$U \oplus V = \{ 0, 1, 3, 8, 9, 11, 12, 13, 15 \}$
$U \ominus V = \{ 0, 1, 3, 4, 5, 7, 12, 13, 15 \}$

今回の実装では、「ずらし」でも、「反射」でもないホモメトリックな解は別々の解として区別している。

アルゴリズム

バイオインフォマティクスのためのアルゴリズム入門には３つのアルゴリズムが載っている。各アルゴリズムは、$(\Delta X, |X|)$ を入力として取る。

全探索アルゴリズムの例1
全探索アルゴリズムの例2
Skiena のアルゴリズム

1. 全探索アルゴリズムの例1

実装: BruteForcePDP.go

大雑把なアルゴリズム。問題の本質とはあまり関係のない ${\max(\Delta X)}$ の値が大きいと組み合わせの数が膨大になり遅い。

$1,2,\dots,\max(\Delta X) - 1$ の中から重複無く $|X| - 2$ 個選んだ時の組み合わせ $x_1,x_2,\cdots,x_{|X|-2}$ を列挙し、列挙した各組み合わせについて、$\Delta \{0,x_1,x_2,\cdots,x_{|X|-2},\max(\Delta X)\} = \Delta X$ であれば、$\{0,x_1,x_2,\cdots,x_{|X|-2},\max(\Delta X)\}$ を解とする。

なお、Go で 0 から始まる連続する n 個の整数を重複無く k 個選んだ時の組み合わせの列挙 - Qiitaを書いていたのはこのアルゴリズムのためだったりもする

2. 全探索アルゴリズムの例1

実装: AnotherBruteForcePDP.go

上記のアルゴリズムを少し改良したもの。$1,2,\dots,\max(\Delta X) - 1$ の中から選ぶのではなく、$\Delta X$ の中から選ぶようする。${\max(\Delta X)}$ に振り回されることはなくなるが、 $|X|$ が大きくなると組み合わせの数が膨大になり、やっぱり遅い。

3. Skiena のアルゴリズム

実装: SkienaPDP.go

元の論文は多分これ

A Partial Digest Approach to Restriction Site Mapping
Steven Skiena
- Reconstructing Sets from Interpoint Distances

下記リンク先のPDFに載ってる擬似コードがわかりやすい。

Past Projects and Papers
- Partial Digest Problem (Bioinformatics)

なお、リンク先は解が一つ見つかればそこで打ち切る擬似コードだが、今回の実装では下記のように区別している解を全て列挙している。

一般性を失うこと無く、求める解の左端を 0, 右端を $width \leftarrow \max(\Delta X)$ と置くことができる。 $X' \leftarrow \{0, width\}$
$\Delta X \setminus \Delta X' = \emptyset$ であれば $X'$ が解
$l_{\max} \leftarrow \max (\Delta X \setminus \Delta X')$ とし、下記二通りの可能性を試す
- 左端から $l_{\max}$ だけ離れた位置にある点（右端より）が解に含まれている可能性
  - $X_r' \leftarrow X' \cup \{l_{\max}\}$ とし、$\Delta X_r' \subset \Delta X$ であれば、$X_r'$ について、2,3 を再帰的に実行
- 右端から $l_{\max}$ だけ離れた位置にある点（左端より）が解に含まれている可能性
  - $X_l' \leftarrow X' \cup \{width - l_{\max}\}$ とし、$\Delta X_l' \subset \Delta X$ であれば、$X_l'$ について、2,3 を再帰的に実行

分岐していくので、最悪の場合（バックトラックが発生しない場合）計算時間が指数関数時間になる。ほとんどの問題のインスタンスは分岐が発生しないので、そこまで悪くはならない。なお、下記論文に指数関数時間になる例が載っているっぽい。

An Exponential Example for a Partial Digest Mapping Algorithm

本に載ってないアルゴリズム

本には「Nivat らのアルゴリズムは多項式時間のPDPアルゴリズムだぜイエーイ」みたいなことが書いてあるが、その脚注には「問題の形式が違うからPDPが多項式時間で解けるかはわからない」といったことが書いてある。謎い。
恐らく、Nivat らのアルゴリズムはThe chords’ problemを指していると思われるが、英文を読む気が起きないので確かめていない。

軽くググッた感じだと、A fast algorithm for the partial digest problemが「爆速だぜ！」感出してる、でも、英文を読む気が起きないので確かめていない。

ベンチマーク

等差数列、等比数列、フィボナッチ数列が解となるテストデータを使ってベンチマークを取ってみた。

等差数列が解になる場合

BenchmarkPDP_ArithmeticSeq_test.go

下記の表はそれぞれ cpu 数 1 と 8 の場合の結果である
{0,2,4,...,18}などは、解となる集合を表す
数値は各実装の１処理当たりの実行時間を表す[ns/op]

-cpu 1	{0,2,4,...,18}	{0,2,4,...,248}	{0,2,4,...,998}
BruteForcePDP	2005732790	×	×
AnotherBruteForcePDP	457334	993430077	×
SkienaPDP	65709	9917569	146312827

-cpu 8	{0,2,4,...,18}	{0,2,4,...,248}	{0,2,4,...,998}
BruteForcePDP	492614671	×	×
AnotherBruteForcePDP	214647	240209341	×
SkienaPDP	113214	27767754	493899901

結果について
- BruteForcePDP が圧倒的に遅く、SkienaPDP が圧倒的に速い
- BruteForcePDP と AnotherBruteForcePDP は cpu 数を増やすと早くなるので、並行化が上手く効いている
- SkienaPDP は cpu 数を増やすと遅くなる。今回の実装では、左右対称な解に対して分岐しないように制御していることが原因だと思われる。

等比数列が解になる場合

BenchmarkPDP_GeometricSeq_test.go

下記の表はそれぞれ cpu 数 1 と 8 の場合の結果である
{0,1,3,7...,2^6-1}などは、解となる集合を表す
数値は各実装の１処理当たりの実行時間を表す[ns/op]

-cpu 1	{0,1,3,7...,2^6-1}	{0,1,3,7...,2^8-1}	{0,1,3,7...,2^60-1}
BruteForcePDP	1041527893	×	×
AnotherBruteForcePDP	9363707	3691641539	×
SkienaPDP	72956	162724	55904758

-cpu 8	{0,1,3,7...,2^6-1}	{0,1,3,7...,2^8-1}	{0,1,3,7...,2^60-1}
BruteForcePDP	296637903	×	×
AnotherBruteForcePDP	2758850	916111740	×
SkienaPDP	60134	119504	25039748

結果について
- BruteForcePDP が圧倒的に遅く、SkienaPDP が圧倒的に速い
- どれも cpu 数を増やすと早くなるので、並行化が上手く効いている
- どれも等比数列が解になる場合よりも、等差数列が解になる場合の方が得意っぽい
  - BruteForcePDP は問題のサイズに対して、等比数列よりも等差数列の方が ${\max(\Delta X)}$ が小さくなるためだと思われる
  - AnotherBruteForcePDP は問題のサイズに対して、等比数列よりも等差数列の方が $\Delta X$ を通常の集合とした時の濃度が小さくなるためだと思われる
  - SkienaPDP は、今回の実装では左右対称な解に対して分岐しないように制御しているためだと思われる。

フィボナッチ数列が解になる場合

BenchmarkPDP_Fib_test.go

下記の表はそれぞれ cpu 数 1 と 8 の場合の結果である
{F_0,F_1,F_3,F_4...,F_8}などは、解となる集合を表す
数値は各実装の１処理当たりの実行時間を表す[ns/op]

-cpu 1	{F_0,F_1,F_3,F_4...,F_8}	{F_0,F_1,F_3,F_4...,F_9}	{F_0,F_1,F_3,F_4...,F_40}
BruteForcePDP	1948687444	×	×
AnotherBruteForcePDP	581294824	15704345634	×
SkienaPDP	172652	228241	16400041

-cpu 8	{F_0,F_1,F_3,F_4...,F_8}	{F_0,F_1,F_3,F_4...,F_9}	{F_0,F_1,F_3,F_4...,F_40}
BruteForcePDP	501265633	×	×
AnotherBruteForcePDP	148893054	3991023144	×
SkienaPDP	114233	152727	7557608

結果について
- SkienaPDP が圧倒的に速い
- どれも cpu 数を増やすと早くなるので、並行化が上手く効いている
- BruteForcePDP と AnotherBruteForcePDP の差が、等比数列が解になる場合や、等差数列が解になる場合ほど顕著でない。フィボナッチ数列が解になる場合 $\Delta X$ を通常の集合とした時の濃度があまり小さくならないためだと思われる

感想

SkienaPDP の並行化について

SkienaPDP が等差数列が解になる場合を除いて、cpu 数を増やすと実効速度がだいたい２倍位になる。これは、左右反転の解を別々の cpu で同時に計算して早くなってる可能性がある。なので、左右反転の解を別々に計算しないように工夫すれば、平行処理がなくても２倍早くなるのではないかという気がしている。

Go おもろい

簡単に並行処理書けるの良い
range便利
- goroutine と channel を Generator 的に使いまくってるけどあまりよろしくないらしい
go test 便利
- seq 1 8 | xargs -n 1 go test -cpu とかで GOMAXPROCS の値を変えながらテストできるの良い
- go test -covermode=count -coverprofile=count.out; go tool cover -func=count.out | grep -v 100 とかで、0% の死んでる関数探せるの良い
- go test -bench hoge で簡単にベンチマーク取れるの良い
- 以下は今後使えたら良いなって思ってる
  - go test -bench . -benchmem
  - go test -covermode=count -coverprofile=count.out; go tool cover -html=cover.out