Rustで順列の生成を行う

概要

巡回セールスマン問題(TSP)で厳密解を求めるために Necklace Permutation を効率よく生成する必要があった。そして、Rustにおいては itertools など順列を生成できるライブラリはあったが、 Necklace Permutation については存在しなかった。

Qiita には Permutation の記事は山ほどあるが、アルゴリズムの設計に役立つ記事はなさそうなので、記事を残しておく。ちなみに、今回は Permutation の記事である。

問題設定

Necklace Permutation を実装するにあたってのウォーミングアップである。

• nPn の形の順列を生成できる
• メモリ使用量が n の定数倍程度におさまること
  • 順列は配列ではなくイテレータで取得できる必要がある
  • 時間がかかっても順列を生成し続けられることが重要
• 順序は規則的でなくともよい (ソートはしなくてよい)
• 置換の符号はいらない

参考: 推定メモリ使用量

イテレータで動作するアルゴリズムは必須ですよね。
TSPLIB は 24～13509 都市まであるが、小規模なものであれば厳密解の分布を調べたいと思っている。でも、15未満ぐらいだろう。

すべての結果を配列に格納する版でのメモリ使用量の影響をまとめる。

メモリ使用量は次の計算で概算する。

(メモリ使用量) = n![個数/経路長] * 4[byte/個] * n [経路長]

n	n!	メモリ使用量(概算)
1	1	4 B
2	2	16 B
3	6	72 B
4	24	384 B
5	120	2 KiB
6	720	16 KiB
7	5040	137 KiB
8	40320	1 MiB
9	362880	12 MiB
10	3628800	138 MiB
11	39916800	1675 MiB
12	479001600	21 GiB
13	6227020800	302 GiB
14	87178291200	4 TiB
15	1307674368000	71 TiB
16	20922789888000	1218 TiB
17	355687428096000	21 PiB
18	6402373705728000	409 PiB
19	121645100408832000	8 EiB
20	2432902008176640000	169 EiB

数学的なおさらい

• 巡回群のように一つの生成元(カウンター)なりで生成できない

任意の置換は互換の積で表すことができる

実際には可換なものもある。
そのため、素朴に置換の組み合わせすべてをプログラムで生成するときは停止しないか、重複のチェックに大容量のメモリ領域が必要で、直接はアルゴリズムにはできない。

\sigma = \sigma_1 \sigma_2  \cdots \sigma_n

$\sigma_1 = (1 ,,,, 2), \sigma_2 = (3 ,,,, 4)$

可換なものもある。

(4 3)(1 2)(1 2 3 4) = (2 1 4 3)
(1 2)(4 3)(1 2 3 4) = (2 1 4 3)

対称群 (Symmetric group)と同型

Symmetric group - Wikipedia より、個数は n! である。また、つぎのような関係など、可換となる演算があるので、生成の仕方は考慮が必要。

• $\sigma_i ^2 = 1$
• $(\sigma_i \sigma_{i+1})^3 = 1$
• など

アルゴリズム検討(ウォーミングアップ)

• 再帰
  • 遅い
  • 規模に応じてメモリ使用量が劇的に増える
  • 遅延評価ができない (イテレータとして使えない)
• 非再帰
  • 再帰より早い
  • メモリ使用量については実装次第
  • 遅延評価ができるかどうかも実装次第

1. 再帰 (引き抜き)

最初に作ってみたもの。遅い。メモリ使用量大きい。

1. [1, 2, 3, 4] などの文字種に対して、[1, 2, 3, 4]から文字を引き抜く
   -> {[2, 3, 4] -> out:[1], [1, 3, 4] -> out:[2], [1, 2, 4] -> out:[3], [1, 2, 3] -> out:[4]}
2. [1, 2, 3, 4] などの文字種に対して、{[1, 2, 3], ...} の文字に対して再帰的に繰り返す
3. 空になったら終了

pub fn gen_perm_with_tmp_vec_u32(v1: Vec<u32>, v2: &mut Vec<u32>, out: &mut Vec<Vec<u32>>)
{
    match v1.len() {
        0 => (),
        1 => {
            v2.push(v1[0]);
            out.push(v2.to_vec());
        },
        _ => {
            for some_x in &v1 {
                let mut vc1 = v1.clone();
                let mut vc2 = v2.clone();
                vc1.retain(|&cur| cur != *some_x);
                vc2.push(*some_x);
                gen_perm_with_tmp_vec_u32(vc1, &mut vc2, out);
            }
        },
    }
}

# [cfg(test)]
mod tests {
    use crate::v1_recursive_ordered_u32::gen_perm_with_tmp_vec_u32;

    #[test]
    fn test_gen_perm_3() {
        let mut result = Vec::<Vec<u32>>::new();
        gen_perm_with_tmp_vec_u32(vec![1, 2, 3], &mut vec![], &mut result);
        println!("{:?}", result);
        assert_eq!(result, vec![
            vec![1, 2, 3],
            vec![1, 3, 2],
            vec![2, 1, 3],
            vec![2, 3, 1],
            vec![3, 1, 2],
            vec![3, 2, 1],
        ]);
    }
}

2. 再帰 (swap)

引き抜きを swap に置き換えたもの。引き抜き版より、少し早い。
メモリ使用量は少し改善しているが依然として大きい。

pub fn gen_perm_with_depth(v: Vec<u32>, m: usize, out: &mut Vec<Vec<u32>>)
{
    if m == v.len() {
        out.push(v.to_vec());
        return;
    }

    for i in m .. v.len() {
        let mut v_new = v.clone();
        if i != m {
            v_new.swap(m, i);
        }
        gen_perm_with_depth(v_new, m+1, out);
    }
}

# [cfg(test)]
mod tests {
    use crate::v2_recursive_unordered_u32::gen_perm_with_depth;

    #[test]
    fn test_gen_perm_3() {
        let mut result = Vec::<Vec<u32>>::new();
        gen_perm_with_depth(vec![1, 2, 3], 0, &mut result);
        println!("{:?}", result);
        assert_eq!(result, vec![
            vec![1, 2, 3],
            vec![1, 3, 2],
            vec![2, 1, 3],
            vec![2, 3, 1],
            vec![3, 2, 1],
            vec![3, 1, 2],
        ]);
    }
}

3. 再帰 (swap; Generics)

Generics版であるが、u32版と比較して速度差はない。

作業用の配列の代わりに再帰の深さを指定している。
逐次版のアイディアの元になっている。

Generics版を作ると where句で型が満たすべき制約を指定することになるのだが、このアルゴリズムでは std::cmp::PartialEq を使っている。同値性がこのアルゴリズムには必要ということである。

pub fn gen_perm_with_depth<T>(v: Vec<T>, m: usize, out: &mut Vec<Vec<T>>)
where T: Clone + std::cmp::PartialEq
{
    if m == v.len() {
        out.push(v.to_vec());
        return;
    }

    for i in m .. v.len() {
        let mut v_new = v.clone();
        if i != m {
            v_new.swap(m, i);
        }
        gen_perm_with_depth(v_new, m+1, out);
    }
}

mod tests {
    use crate::v3_recursive_unordered_gen::gen_perm_with_depth;

    #[test]
    fn test_gen_perm_3() {
        let mut result = Vec::<Vec<u32>>::new();
        gen_perm_with_depth(vec![1, 2, 3], 0, &mut result);
        println!("{:?}", result);
        assert_eq!(result, vec![
            vec![1, 2, 3],
            vec![1, 3, 2],
            vec![2, 1, 3],
            vec![2, 3, 1],
            vec![3, 2, 1],
            vec![3, 1, 2],
        ]);
    }
}

4. 逐次版 (swap)

再帰版より、倍近く早い。アイディアは再帰版に近い。
まだ配列で持っているので、メモリ使用量は依然として大きい。
result.push の部分で結果を覚えておく必要があり、出力用の配列がまだ必要。このアルゴリズムは、重複がないことが前提となっている。なお、 std::cmp::PartialEq を使っていないが比較するロジックも加えれば、重複の場合も対応できるだろう。

1. 1文字目のバリエーションを洗い出し尽くして、結果リストに登録
2. 結果リストから取り出したものに対して、2文字目のバリエーションを洗い出してリストに登録
3. 上記を最後の文字まで続ける

pub fn gen_perm<T>(v: Vec<T>)
-> Vec<Vec<T>>
where T: Clone
{
    let num_of_chars = v.len();
    let mut result = Vec::<Vec<T>>::new();
    result.push(v);
    for n in 0 .. num_of_chars {
        let result_len = result.len();
        for result_idx in 0..(result_len) {
            for i in (n+1) .. num_of_chars {
                let mut v_new = result[result_idx].clone();
                v_new.swap(n, i);
                result.push(v_new);
            }
        }
    }
    result
}

# [cfg(test)]
mod tests {
    use crate::v4_iterative_unordered::gen_perm;

    #[test]
    fn test_gen_perm_3() {
        let result = gen_perm(vec![1, 2, 3]);
        println!("{:?}", result);
        assert_eq!(result, vec![
            vec![1, 2, 3],
            vec![2, 1, 3],
            vec![3, 2, 1],
            vec![1, 3, 2],
            vec![2, 3, 1],
            vec![3, 1, 2],
        ]);
    }
}

5. 逐次版 (swap, ソートあり)

結果が正しいのか調べる趣旨でソート版も作った。
処理は 4. とほぼ同じで処理が終わった後に sort_by で並び替えている。
結果は、再帰版の1.5倍ぐらい遅い。

pub fn gen_perm<T>(v: Vec<T>)
-> Vec<Vec<T>>
where T: Clone + std::cmp::PartialOrd
{
    let num_of_chars = v.len();
    let mut result = Vec::<Vec<T>>::new();
    result.push(v);
    for n in 0 .. num_of_chars {
        let result_len = result.len();
        for result_idx in 0..(result_len) {
            for i in (n+1) .. num_of_chars {
                let mut v_new = result[result_idx].clone();
                v_new.swap(n, i);
                result.push(v_new);
            }
        }
    }
    result.sort_by(|a, b| {
        let m = a.len() - 1;
        for i in 0 .. m {
            if a[i] != b[i] {
                return a[i].partial_cmp(&b[i]).unwrap();
            }
        }
        return a[m].partial_cmp(&b[m]).unwrap();
    });
    result
}

アルゴリズム実装 (イテレータ版)

解説

「4. 逐次版 (swap)」より、各階層の文字種の数だけループすることで、期待する結果が得られることが分かっている。期待する結果とは、全パターンを網羅でき、抽出済みのものが出てこない置換が行えることである。

例えば、n=4 であれば、ループ回数は次のようになる。文字数が減るごとにループ回数も減る。

n文字目	loop回数
1	4
2	3
3	2
4	1

上記より、4文字目のループは1回なので不要ではあるが、4次元の整数ベクトルで組み合わせを一意に表現できることを意味する。よって、イテレータとしては、上記のループ回数と生成前の入力を覚えておく領域があればよい。

コード

• new の部分で、n文字種のループ用の変数とループ範囲の配列を作っている
  • ループ回数については、毎ループごとにカウンターが最大になったかどうかを判定することが無駄が大きく、桁あふれ用の領域をもつ(11!で4000万回動作するので if 一つも数%の影響がある)
• next ではカウンターの加算処理とカウンター溢れの終了の制御を行っている
• 毎度カウンターに該当する置換をしているが、キャッシュする処理があればさらなる高速化ができるかもしれない

use std::ops::Range;

struct PermutationIterator<T> {
    initial: Vec<T>,
    ranges: Vec<Range<u16>>,
    indexes: Vec<u16>,
}

impl<T> PermutationIterator<T> {
    #![allow(dead_code)]
    pub fn new(p: Vec<T>) -> PermutationIterator<T>
    where T: Clone
    {
        let mut indexes = Vec::<u16>::with_capacity(p.len());
        let mut ranges = Vec::<Range<u16>>::with_capacity(p.len());
        indexes.push(0 as u16);
        ranges.push(0 as u16 .. 1 as u16);
        for i in 1 .. p.len() {
            ranges.push((i-1) as u16 .. p.len() as u16);
            indexes.push(ranges[i].start.clone());
        }
        PermutationIterator {
            initial: p.clone(),
            ranges: ranges,
            indexes: indexes,
        }
    }
}

impl<T: Clone> Iterator for PermutationIterator<T> {
    type Item = Vec<T>;
    fn next(&mut self) -> Option<Self::Item>
    {
        if self.indexes[0] > 0 {
            return None;
        }
        let mut v = self.initial.clone();
        let end = self.initial.len();
        for i in 0 .. (end-1) {
            if i != (self.indexes[i+1] as usize) {
                v.swap(i, self.indexes[i+1] as usize);
            }
        }
        self.indexes[end - 1] += 1;
        for i in (1 .. end).rev() {
            if self.indexes[i] >= self.ranges[i].end {
                self.indexes[i] = self.ranges[i].start;
                self.indexes[i - 1] += 1;
            }
        }
        return Some(v);
    }
}

# [cfg(test)]
mod tests {
    #[test]
    fn test_perm_iter_3() {
        use crate::PermutationIterator;

        let mut iter = PermutationIterator::new(
            vec![1 as u8, 2, 3]
        );
        assert_eq!(iter.next(), Some(vec![1, 2, 3]));
        assert_eq!(iter.next(), Some(vec![1, 3, 2]));
        assert_eq!(iter.next(), Some(vec![2, 1, 3]));
        assert_eq!(iter.next(), Some(vec![2, 3, 1]));
        assert_eq!(iter.next(), Some(vec![3, 2, 1]));
        assert_eq!(iter.next(), Some(vec![3, 1, 2]));
        assert_eq!(iter.next(), None);
    }
}

まとめ

• 再帰をベースとするアイディアはメモリ使用量が大きくなる傾向がある
• イテレータでも順列の生成ができる
  Rustの itertools も似たようなことをしているが、状態管理がすごく複雑で読み解くことが難しい)
• nが大きくなってくると、重複あり/なしの前提も速度に影響がでてくる
  • 重複なしなら、要素同士の等値性や比較なしで生成することもできる
• 何次元で生成された値が表現できるだろうか…と考えるとイテレータ版のアプローチが見つかることもある (再帰版だとあまりそういう考えにならない)
• Rust の Generics は数学的に面白い
  • 型に求められる最低限な演算を意識させてくれる

関連記事

• 2つの操作のみで全順列を列挙する：対称群のグラフ上のハミルトン路にもとづく順列生成の紹介と実装 : サルノオボエガキ
  これと似たようなことを紙上で試していた。
  サイクルがあるとメモリ効率の良いアルゴリズムを作るのが難しい。
• 再帰を使わない順列生成 : 1から9の数字を1回ずつ使ってできる数を、小さい方から順に列挙する(アルゴリズム編) - Javaエンジニア、React+Firebaseでアプリを作る
• rust-itertools / itertools