buyoh.hateblo.jp

残念な競技プログラミング参加者による学習記録

最小頂点被覆の指数時間アルゴリズムの実装(その1,半分全列挙)

2018/02/22 軽く更新.

その1.

最終目標

これをACする.

  • CODE THANKS FESTIVAL 2017 - G Mixture Drug

code-thanks-festival-2017-open.contest.atcoder.jp

知識

最小頂点被覆問題ってなんぞ?

グラフG=(V,E)が与えられる.
頂点集合S⊆Vを考える.
全ての辺e∈Eが,Sに含まれる頂点に<u>1つ以上</u>接続している時,SはグラフGの頂点被覆と呼ぶ.

集合Sの大きさが最小であるようなグラフGの頂点被覆Sを求め,その大きさを出力してください.

note

  • 頂点数N,辺の数M
  • 与えられたグラフは単純グラフである
  • 辺は1本以上存在する
  • 頂点に割り当てられた重みは全て等しい(考慮しない).

と仮定します.

説明で使うグラフ構造体を示します.

struct Graph {
    size_t n;
    vector<vector<int>> vertex_to;

    Graph(size_t n = 1) :n(n), vertex_to(n) {}
    void connect(int from, int to) {
        vertex_to[(size_t)from].emplace_back(to);
        vertex_to[(size_t)to].emplace_back(from);
    }
};

基礎

  • 工夫する必要が無いのであれば,時間計算量 O((N+M)2N),空間計算量 O(N) の次のような全探索が思いつく.
  • choice[i]true の時,その頂点が選ばれることを示す.
int vertex_cover(const Graph& graph) {
    const int n = graph.n;

    int best_bit = 0;
    int best_cnt = 1e9;

    repeat(bit, 1<<n) {
        bitset<31> choice = bit;
        bool ok = true;
        repeat(j, n) {
            if (choice[j] == false) {
                for (int to : graph.vertex_to[j])
                    if (!choice[to]) {
                        ok = false; j = n;
                        break;
                    }
            }
        }
        // bitが頂点被覆ならば
        if (ok) {
            int cnt = choice.count();
            if (cnt < best_cnt) {
                best_cnt = cnt;
                best_bit = bit;
            }
        }
    }
    return best_cnt;
}

動的計画法 + 半分全列挙による高速化

準備(動的計画法(通称:bitDP))

  • 半分全列挙の説明をする前に,次のような問題を考えたい.
グラフG=(V,E)が与えられる.
全てのS⊆Vについて,GからV-Sを取り除いてできたグラフG'の最小頂点被覆のサイズを求めよ.
  • この問題は時間計算量 O( (N+M)2N) の動的計画法 *1 で求まる.
  • dp[S] を S⊆Vの頂点被覆のサイズとする.
  • dp[S] が独立集合ならば,グラフには辺は含まれていない.
    • 頂点被覆のサイズdp[S] = 0 とすることができる.
  • Sに頂点v∉Sを加えたS'を考える.
    • dp[S']は,dp[S]+1を超えないはず.
    • dp[S'] = min(dp[S'], dp[S] + 1) とすることができる.
  • 以上から,次のように実装できる.
int vertex_cover(const Graph& graph) {
    const int n = graph.n;
    const int n_half = n/2;
    const int inf = 1e9;

    // dp[S] : S \subseteq V, vertexcover of G[S] 
    vector<int> ans_partial(1 << n, inf);
    ans_partial[0] = 0;

    for (int bit = 1; bit < 1 << n; ++bit) {
        // S
        bitset<31> mask = bit;

        bool independent = true;
        for (int i = 0; i < n; ++i) {
            if (mask[i] == false) continue;
            for (int j : graph.vertex_to[i])
                if (mask[j] == true) {
                    independent = false;
                    i = n; break;
                }
        }
        // maskが独立集合なら,最小頂点被覆は0
        if (independent == true) {
            ans_partial[bit] = 0;
        }
        // 遷移
        for (int i = 0; i < n; ++i) {
            if (mask[i] == true) continue;
            ans_partial[bit | (1 << i)] = min(
                ans_partial[bit | (1 << i)],
                ans_partial[bit] + 1
            );
        }
    }

    return ans_partial[(1 << n) - 1];
}

半分全列挙

  • 半分全列挙を思い出そう.
    1. 集合を A,B の2つに分ける.
    2. Aの部分集合の解をすべて求めておく.
    3. Bのすべての部分集合を列挙.S⊆Bについて,適切なAの部分集合の解を組み合わせる.
    4. 最もナイスな解を出力. .
  • 最小頂点被覆における半分全列挙を考えてみる.
    1. 頂点集合 V を 集合 A,B に適当に分ける.
    2. Aのすべての部分集合について,最小頂点被覆を求める.
    3. Bのすべての部分集合を列挙.
      1. S⊆B が頂点被覆でないなら解ではない.
      2. 適切なAの部分集合の解を組み合わせる.
    4. 最もナイスな解を出力. .
  • ところで,3.2. の適切なAの部分集合とは?
    • 例えば,何も考えず適当にAの最小頂点被覆の解を取り出したらどうなるか.
      • Wrong Answer.なぜならば,A,B のカットエッジ *2 の被覆をチェックしていないから.
      • 逆に,カットエッジ以外の辺は 1. と 3.1 で被覆されていることが保証される.
    • B によって被覆されていないカットエッジの要素を求め,そのような辺があれば A 側で必ず選択するようにする.
int vertex_cover(const Graph& graph) {
    const int n = graph.n;
    const int n_A = n / 2;
    const int n_B = n - n_A;
    const int inf = 1e9;

    // 少なすぎるケースは回避
    if (n < 6) return vertex_cover_l(graph);

    // group A : i <  n_A
    // group B : i >= n_A

    // dp[S] : S \subseteq V_A, 
    vector<int> ans_partial_A(1 << n_A, inf);
    ans_partial_A[0] = 0;

    // groupA の部分集合の解を全列挙
    // ここでは groupA と groupB のカットエッジの被覆は考えない.
    for (int bit = 1; bit < 1 << n_A; ++bit) {
        bitset<31> mask = bit;

        bool independent = true;
        for (int i = 0; i < n_A; ++i) {
            if (mask[i] == false) continue;
            for (int j : graph.vertex_to[i])
                if (j < n_A && mask[j] == true) {
                    independent = false;
                    i = n_A; break;
                }
        }
        // maskが独立集合なら,最小頂点被覆は0
        if (independent == true) {
            ans_partial_A[bit] = 0;
        }
        // 遷移
        for (int i = 0; i < n_A; ++i) {
            if (mask[i] == true) continue;
            ans_partial_A[bit | (1 << i)] = min(
                ans_partial_A[bit | (1 << i)],
                ans_partial_A[bit] + 1
            );
        }
    }

    // answer
    int best = inf;

    // groupB の部分集合を全列挙
    for (int bit = 0; bit < 1 << n_B; ++bit) {
        bitset<31> choice_B = bit;

        bool cover = true;
        for (int _i = 0; _i < n_B; ++_i) {
            if (choice_B[_i] == true) continue;
            for (int j : graph.vertex_to[_i + n_A])
                if (n_A <= j && choice_B[j - n_A] == false) {
                    cover = false;
                    _i = n_A; break;
                }
        }
        // choice_Bが頂点被覆でないなら,reject
        if (cover == false) continue;

        bitset<31> mask_A = (1 << n_A) - 1;
        int fix_A = 0;

        for (int _i = 0; _i < n_B; ++_i) {
            if (choice_B[_i] == true) continue;
            for (int j : graph.vertex_to[_i + n_A]) {
                // choice_B で選んでいない頂点が,groupAと隣接するならば
                if (j < n_A && mask_A[j]) {
                    // そのgroupAの頂点は必ず選択する.
                    mask_A[j] = 0; // 頂点jは考慮しなくてよい
                    ++fix_A;       // 必ず選択することにしたので
                }
            }
        }

        best = min(best, ans_partial_A[mask_A.to_ullong()] + (int)choice_B.count() + fix_A);
    }

    return best;
}
  • 頂点を半分だけ全列挙するのに掛かる時間は O(2N/2)
  • 独立集合・頂点被覆の判定に掛かる時間は O(N+M)
  • カットエッジの本数は最大でも O(M)
  • 以上により,時間計算量はO((N+M)2N/2)
  • 半分全列挙したものを保持する必要があるので,空間計算量は O(2N/2)
  • Submission #1854558 - CODE THANKS FESTIVAL 2017(Parallel) | AtCoder

  • 時間計算量の M は取り除けないだろうか?

    • A の半分全列挙時の独立集合判定と,カットエッジの被覆判定でO(M)がでてくる.

独立集合判定のO(M)の除去

  • あらかじめ全ての頂点部分集合について,独立集合判定をしよう.
  • 頂点数が0の時,yes.
  • 頂点が1つだけの時,yes.
  • 頂点が2つの時,間に辺があれば no.
  • 3つ以上の時は漸化式で.
    • S⊆V, T⊆V が既に分かっていて,S∪T が分からないとき,
    • DP[S∪T] = DP[S] 論理積 DP[T]
  • 時間 O(M+N2N/2) で計算可能.
int vertex_cover(const Graph& graph) {
    /* 略 */

    // dp[S] : S \subseteq V_A, 
    vector<int> ans_partial_A(1 << n_A, 0);

    // Sが独立集合なら0,独立集合でないならinfとなるように初期状態を構成する.
    for (int i = 0; i < n_A; ++i) {
        int bit = 1 << i;
        for (int j : graph.vertex_to[i])
            if (j < n_A)
                ans_partial_A[bit | (1 << j)] = inf; // 頂点数2
    }
    for (int bit = 3; bit < 1 << n_A; ++bit) {
        if (bitset<31>(bit).count() <= 1) continue;
        for (int i = 0; i < n_A; ++i) {
            ans_partial_A[bit | (1 << i)] |= ans_partial_A[bit];
        }
    }

    // groupA の部分集合の解を全列挙
    // ここでは groupA と groupB のカットエッジの被覆は考えない.
    for (int bit = 1; bit < 1 << n_A; ++bit) {
        bitset<31> mask = bit;

        // (略)
        // // maskが独立集合なら,最小頂点被覆は0
        // if (independent == true) {
        //     ans_partial_A[bit] = 0;
        // }
        // 遷移
        for (int i = 0; i < n_A; ++i) {
            if (mask[i] == true) continue;
            ans_partial_A[bit | (1 << i)] = min(
                ans_partial_A[bit | (1 << i)],
                ans_partial_A[bit] + 1
            );
        }
    }
    /* 略 */
}
  • 独立集合の判定を追い出すことができた.
  • 上の説明とは異なる実装をしている.
    • yes=0,no=inf としたので,論理積論理和に置き換わっている.
    • このような実装にすることで,そのまま半分全列挙のDPに流用できる.

独立集合判定のO(M)の除去

  • O(M) はどこで出現しているか?
    • choice_B で選んでいない頂点と隣接する,groupA の頂点の列挙に O(M) 掛かる
  • 漸化式・動的計画法を用いて一括で列挙してみよう.

  • つまり,次のようなDP配列を求めれば良さそう.

全ての β⊆B について,βに隣接する頂点の集合 α⊆A を求める.
  • βの頂点数が0,1の時,簡単に求まる(隣接頂点を列挙).
  • 2つ以上の時は漸化式で.
    • S⊆B, T⊆B が既に分かっていて,S∪T が分からないとき,
    • DP[S∪T] = DP[S] ∪ DP[T]
    • 和集合の演算時間は O(N).
  • よって,O(N2N) で全部求まる.
    // dp[S] Sに隣接する頂点のbit配列
    vector<int> adjacent2B(1 << n_B, 0);
    for (int _i = 0; _i < n_B; ++_i) {
        bitset<31> bits = 0;
        for (int j : graph.vertex_to[n_A + _i])
            if (j < n_A)
                bits[j] = 1;
        adjacent2B[1 << _i] = (int)bits.to_ulong();
    }
    for (int bit = 1; bit < 1 << n_B; ++bit) {
        for (int _i = 0; _i < n_B; ++_i) {
            adjacent2B[bit | (1 << _i)] |= adjacent2B[bit];
        }
    }
  • これを組み込む.
int vertex_cover(const Graph& graph) {

    /* 略 */

    // dp[S] Sは独立集合ではない
    vector<int> not_independent_B(1 << n_B, 0);

    // dp[S] Sに隣接する頂点のbit配列
    vector<int> adjacent2B(1 << n_B, 0);

    /* 略 */

    // answer
    int best = inf;

    // groupB の部分集合を全列挙
    for (int bit = 0; bit < 1 << n_B; ++bit) {
        bitset<31> choice_B = bit;

        // choice_Bが頂点被覆でないなら,reject
        if (not_independent_B[((1 << n_B) - 1) ^ bit]) continue;

        bitset<31> mask_A = (1 << n_A) - 1;
        int fix_A = 0;

        bitset<31> adj = adjacent2B[((1 << n_B) - 1) ^ bit];
        for (int i = 0; i < n_A; ++i) {
            // choice_B で選んでいない頂点が,groupAと隣接するならば
            if (adj[i] == true) {
                // そのgroupAの頂点は必ず選択する.
                mask_A[i] = 0; // 頂点jは考慮しなくてよい
                ++fix_A;       // 必ず選択することにしたので
            }
        }
        best = min(best, ans_partial_A[mask_A.to_ullong()] + (int)choice_B.count() + fix_A);
    }

    return best;
}
  • 以上により,時間計算量 O(N2N/2 ) を達成.
  • 定数倍がクソ重いので 改善後 の実行時間が 改善前 の半分にもなっていないのが面白いです

成果物

inline int bitcount(int x) {
    return bitset<31>(x).count();
}

// 半分全列挙
int vertex_cover(const Graph& graph) {
    const int n = graph.n;
    const int n_A = n / 2;
    const int n_B = n - n_A;
    const int inf = 1e9;

    // 少なすぎるケースは回避
    if (n < 6) return vertex_cover_l(graph);

    // group A : i <  n_A
    // group B : i >= n_A

    // dp[S] : S \subseteq V_A, 
    vector<int> ans_partial_A(1 << n_A, 0);

    // Sが独立集合なら0,独立集合でないならinfとなるように初期状態を構成する.
    for (int i = 0; i < n_A; ++i) {
        int bit = 1 << i;
        for (int j : graph.vertex_to[i])
            if (j < n_A)
                ans_partial_A[bit | (1 << j)] = inf; // 頂点数2
    }
    for (int bit = 3; bit < 1 << n_A; ++bit) {
        if (bitcount(bit) <= 1) continue;
        for (int i = 0; i < n_A; ++i) {
            ans_partial_A[bit | (1 << i)] |= ans_partial_A[bit];
        }
    }

    // groupA の部分集合の解を全列挙
    // ここでは groupA と groupB のカットエッジの被覆は考えない.
    for (int bit = 1; bit < 1 << n_A; ++bit) {
        for (int i = 0; i < n_A; ++i) {
            if (bit & (1 << i)) continue;
            ans_partial_A[bit | (1 << i)] = min(
                ans_partial_A[bit | (1 << i)],
                ans_partial_A[bit] + 1
            );
        }
    }

    // dp[S] Sは独立集合ではない
    vector<int> not_independent_B(1 << n_B, 0);
    for (int _i = 0; _i < n_B; ++_i) {
        int bit = 1 << _i;
        for (int j : graph.vertex_to[n_A + _i])
            if (n_A <= j)
                not_independent_B[bit | (1 << (j - n_A))] = 1;
    }
    for (int bit = 3; bit < 1 << n_B; ++bit) {
        if (bitcount(bit) <= 1) continue;
        for (int _i = 0; _i < n_B; ++_i) {
            not_independent_B[bit | (1 << _i)] |= not_independent_B[bit];
        }
    }

    // dp[S] Sに隣接する頂点のbit配列
    vector<int> adjacent2B(1 << n_B, 0);
    for (int _i = 0; _i < n_B; ++_i) {
        int bit = 0;
        for (int j : graph.vertex_to[n_A + _i])
            if (j < n_A)
                bit |= 1 << j;
        adjacent2B[1 << _i] = bit;
    }
    for (int bit = 1; bit < 1 << n_B; ++bit) {
        for (int _i = 0; _i < n_B; ++_i) {
            adjacent2B[bit | (1 << _i)] |= adjacent2B[bit];
        }
    }

    // answer
    int best = inf;

    // groupB の部分集合を全列挙
    for (int bit = 0; bit < 1 << n_B; ++bit) {

        // choice_Bが頂点被覆でないなら,reject
        if (not_independent_B[((1 << n_B) - 1) ^ bit]) continue;

        int mask_A = (1 << n_A) - 1;
        int fix_A = 0;

        int adj = adjacent2B[((1 << n_B) - 1) ^ bit];
        for (int i = 0; i < n_A; ++i) {
            // choice_B で選んでいない頂点が,groupAと隣接するならば
            if (adj & (1 << i)) {
                // そのgroupAの頂点は必ず選択する.
                mask_A ^= 1 << i; // 頂点jは考慮しなくてよい
                ++fix_A;          // 必ず選択することにしたので
            }
        }

        best = min(best, ans_partial_A[mask_A] + bitcount(bit) + fix_A);
    }

    return best;
}

感想

  • __builtin_popcount は置いといて, bitset<31>(x).count() も結構高速なんですね

参考サイト

*1:bitDP,配るDP.

*2: 一方の端点が A に属する頂点,もう一方が B であるような辺の集合