ビットDP(bit DP)の考え方 ~集合に対する動的計画法~
競技プログラミングで良く使われる動的計画法の1種、「ビットDP」と呼ばれるものについてまとめました。
ビットDPとは
ビットDP(bit DP) とは、ビットで表現した集合を添え字に持つ動的計画法(DP)のことです。
基本的には、以下のような DPを考えます。
\(dp[ S ]\) := 部分集合 S に対して\(|S|!\) 通りの順序の中から最適なものを選んだときの、何かしらの値
漸化式の更新式としては、
$$dp[ S \cup \{v\} ] = dp[ S ] + cost(S, v)$$
のように集合を1つずつ増やしていく形になることが多いです。
この集合に対するDPによって、\(n\) 個のものに対して\(n!\) 通りの順序の中から最適なものを計算するのを高速化できる場合があります。
集合をビットで表現する
実際には集合を配列の引数にするために、集合を2進数で表現します。全体集合 {0,1,2,…,n-1} の部分集合 S を n 桁の2進数で表現するのです。
例えば、全体集合 {0,1,2,…,9} の部分集合 {0,3,5} については、0bit目と3bit目と5bit目が 1で、それ以外が0となるような2進数の数字 \(10010100000_2\) で表現することができます。
このような「集合をビットで管理する」という考え方が、ビットDPと呼ばれる理由です。
巡回セールスマン問題
頂点数 n の重み付き有向グラフが与えられる。以下を満たす最短経路の距離を求めよ。
- 閉路である
- 始点と終点を除いて、全ての頂点を1度ずつ通る
リンク：AOJ Traveling Salesman Problem
ビットDPで効率的に計算できる代表例として、以上のような巡回セールスマン問題があります。始点と終点を除いて、全ての頂点を1度ずつ通る閉路のことを、ハミルトン閉路などと言います。
先程の概要だけではイマイチ理解しにくいと思います。実際にこの問題を通して、ビットDPの使い方を見ていきましょう。
全探索では間に合わない場合がある
純粋に全探索することを考えると、条件を満たすような経路は \(n!\) 通り考えることができます。
これを全探索して最短経路を求めても良いのですが、\(n\)が非常に小さくとも、\(n!\) は非常に大きい値になり計算が間に合わないことがあります。
ビットDPで高速化
そこでビットDPを使って高速化することを考えてみましょう。ビットDPの基本的な式から考えると、
- \(dp[ S ]\) := {0,1,2,…,n-1} の部分集合 S を巡回する \(|S|!\) 通りの経路のうち最短のものの距離
としたくなりますが、これでは動的計画法の漸化式を立てることができません。漸化式を立てるためには、「直前にどの頂点にいたのか」という情報が必要不可欠だからです。
そこで以下のように DP を定義して更新式を考えましょう。
\(dp[ S ][ v ]\) := {0,1,2,…,n-1} の部分集合 S を巡回する \(|S|!\) 通りの経路のうち最短のものの距離。ただし、最後に頂点 v に到達した時のみを考える。
更新式は、頂点 u から v までの距離を cost(u,v) とすると以下の通りとなります。
$$dp[ S \cup \{v\} ][ v ] = \min_{u\in S}(dp[ S ][ u ] + cost(u, v))$$
計算量の変化
このようにすると、全探索のときは \(O(n!)\) だった計算量が、 \(O(n^2 2^n)\) 程度にまでなります。
\(O(n!)\) では \(n=10\) を超えるとかなり時間がかかりますが、\(O(n^2 2^n)\)なら軽い計算で \(n=20\) 前後まで数秒で計算ができるようになります。
実装上の注意
最短経路を求める問題なので、存在しないパスは非常に大きい値(INFなど)にしてしまえば大丈夫です。（オーバーフローには注意しましょう。）
また、更新の順序に関しては「集合の要素が1の時, 2の時, 3の時, … 」のように考えるとループで表現することが難しいです。このような時はメモ化再帰などで実装するのが一般的です。
しかし、実際には2進数の数 \(i, j\) が表す集合が \(S(i) \subset S(j)\) ならば \(i \leq j\) となる性質から、以下のようなループで更新してしまっても問題ありません。
for(int S=0; S<(1<<n); S++){
// 更新処理
}
さらに、どの頂点からスタートしても最短となるハミルトン閉路が存在することを考えると、自分で適当にスタートとなる頂点を決めてしまって大丈夫です。
始点と終点を一致させるためには、頂点0からスタートしたとして以下のようにすれば良いです。
dpの定義
- \(dp[ S ][ v ]\) := 頂点0からスタートして、{0,1,2,…,n-1} の部分集合 S を巡回する \(|S|!\) 通りの経路のうち最短のものの距離。ただし、最後に頂点 v に到達した時のみを考える。
dpの初期化
- \(dp[{0}][0] = 0\)
- その他は全て INF
頂点0 からスタートするので以上の通りです。
dpの更新
- \(dp[ S \cup {v} ][ v ] = \min_{u\in S}(dp[ S ][ u ] + cost(u, v))\)
更新式は先程説明した通りです。
答え
- \(dp[\{0,1,…,n-1\}][ 0 ]\)
最終的に頂点0 に到達した時が求める答えになります。
実装例
メモ化再帰での実装と、配るDPを利用した実装です。
メモ化再帰
実装上の注意でも述べましたが、更新の順序に関しては「集合の要素が1の時, 2の時, 3の時, … 」のように考えるとループで表現することが難しいです。
このような時は以下のようにメモ化再帰などで実装するのが一般的です。
#include <bits/stdc++.h>
#define LOOP(n) for (int _i = 0; _i < (n); _i++)
#define REP(i, n) for (int i = 0; i < (n); ++i)
#define RREP(i, n) for (int i = (n); i >= 0; --i)
#define FOR(i, r, n) for (int i = (r); i < (n); ++i)
#define ALL(obj) begin(obj), end(obj)
using namespace std;
using ll = long long;
const int INF = 1000100100;
template <class T>
bool chmin(T &a, const T &b) {
if (b < a) {
a = b;
return 1;
}
return 0;
}
int V, E;
int G[20][20]; // グラフ
int dp[50000][20];
// メモ化再帰
int rec(int S, int v) {
if (S == 0) {
if (v == 0) {
return 0;
} else {
return INF;
}
}
if ((S & (1 << v)) == 0) { // Sに{v}が含まれていない
return INF;
}
int &ret = dp[S][v];
if (ret != 0) return ret;
ret = INF;
REP(u, V) { chmin(ret, rec(S ^ (1 << v), u) + G[u][v]); }
return ret;
}
int main() {
cin >> V >> E;
// グラフの初期化
REP(i, 20) {
REP(j, 20) { G[i][j] = INF; }
}
REP(i, E) {
int s, t, d;
cin >> s >> t >> d;
G[s][t] = d;
}
int ans = rec((1 << V) - 1, 0);
if (ans != INF) {
cout << ans << endl;
} else {
cout << -1 << endl;
}
return 0;
}
配るDP
実際には2進数の数 \(i, j\) が表す集合が \(S(i) \subset S(j)\) ならば \(i \leq j\) となる性質から、配るDPのループを利用した更新順序でも問題なく更新ができます。
#include <bits/stdc++.h>
#define LOOP(n) for (int _i = 0; _i < (n); _i++)
#define REP(i, n) for (int i = 0; i < (n); ++i)
#define RREP(i, n) for (int i = (n); i >= 0; --i)
#define FOR(i, r, n) for (int i = (r); i < (n); ++i)
#define ALL(obj) begin(obj), end(obj)
using namespace std;
using ll = long long;
const int INF = 1000100100;
template <class T>
bool chmin(T &a, const T &b) {
if (b < a) {
a = b;
return 1;
}
return 0;
}
int V, E;
int G[20][20]; // グラフ
int dp[50000][20];
int main() {
cin >> V >> E;
// グラフの初期化
REP(i, 20) {
REP(j, 20) { G[i][j] = INF; }
}
REP(i, E) {
int s, t, d;
cin >> s >> t >> d;
G[s][t] = d;
}
// dp の初期化
REP(i, 50000) {
REP(j, 20) { dp[i][j] = INF; }
}
dp[0][0] = 0; // スタートを頂点 0 とする
REP(S, 1 << V) {
REP(v, V) {
REP(u, V) {
if (S != 0 && !(S & (1 << u))) continue; // u を含まない場合を除く
if ((S & (1 << v)) == 0) {
if (v != u) chmin(dp[S | (1 << v)][v], dp[S][u] + G[u][v]);
}
}
}
}
if (dp[(1 << V) - 1][0] != INF) {
cout << dp[(1 << V) - 1][0] << endl;
} else {
cout << -1 << endl;
}
return 0;
}
G – Revenge of Traveling Salesman Problem
先程のトラベリングセールスマン問題を、重み付き無向グラフについて考える。ただし、パス \(i\) は \(time_i\) 時間経過すると通れなくなる。
この時の最短経路の距離と、最短距離で戻る経路の総数を出力せよ。
リンク：[AtCoder] square869120Contest #1 G – Revenge of Traveling Salesman Problem
先程のトラベリングセールスマン問題の上位互換の問題です。
トラベリングセールスマン問題との相違点
この問題も、基本的にはビットDPによって解くことができますが、いくつか違いがあるので注意しましょう。
有向グラフではなく無向グラフ
ちょっとした違いですが、先程の問題とは異なり無向グラフなので双方向にパスができるように注意しましょう。
入力が非常に大きい
入力が非常に大きく 32bit の整数型では収まらないことに注意しましょう。
最短距離だけではなく、最短距離となる経路の数も求める
最も大きな違いが、最短距離だけではなく、最短距離の経路数まで求める必要があることです。これは以下のように経路数を含めて dp を定義してやればうまく計算ができます。
\(dp[ S ][ v ]\) := {0,1,2,…,n-1} の部分集合 S を巡回する \(|S|!\) 通りの経路のうち最短のものの距離と経路数のペア。ただし、最後に頂点 v に到達した時のみを考える。
パスに時間制限がある
2つ目に大きな違いが、時間によって通れる辺と通れない辺が生じることです。これは更新式を少し変更してやれば対応可能で、
- \(dp[ S \cup \{v\} ][ v ] = \min_{u\in S}(dp[ S ][ u ] + cost(u, v))\) (ただし、\(dp[ S ][ u ] + cost(u, v)\) の値が \(time_{(u,v)}\)以内の時に限る )
などとしてやれば良いです。
もちろん、dp は距離と経路数のペアなので、実際にはもう少し複雑なコードになります。
実装例
メモ化再帰での実装です。
#include <bits/stdc++.h>
#define LOOP(n) for (int _i = 0; _i < (n); _i++)
#define REP(i, n) for (int i = 0; i < (n); ++i)
#define RREP(i, n) for (int i = (n); i >= 0; --i)
#define FOR(i, r, n) for (int i = (r); i < (n); ++i)
#define ALL(obj) begin(obj), end(obj)
using namespace std;
using ll = long long;
const ll INF = 1e16;
struct Edge {
ll d;
ll time;
};
Edge G[16][16]; // グラフ
ll N, M;
pair<ll, ll> dp[1 << 16][16];
pair<ll, ll> rec(ll S, ll v) { // メモ化再帰
if (S == 0) {
if (v == 0) {
return make_pair(0, 1);
} else {
return make_pair(INF, 0);
}
}
if ((S & (1 << v)) == 0) { // Sに{v}が含まれていない
return make_pair(INF, 0);
}
pair<ll, ll> &ret = dp[S][v];
if (ret.first != 0) return ret;
ret = make_pair(INF, 0);
REP(u, N) {
pair<ll, ll> p = rec(S ^ (1 << v), u);
p.first += G[u][v].d;
if (p.first <= G[u][v].time) { // 制限時間以内に通れる時
if (ret.first > p.first) {
ret = p;
} else if (ret.first == p.first) { // 距離が同じ時は経路数を合わせる
ret.second += p.second;
}
}
}
return ret;
}
int main() {
cin >> N >> M;
// グラフの初期化
REP(i, 16) {
REP(j, 16) {
G[i][j].d = INF;
G[i][j].time = 0;
}
}
REP(i, M) {
ll s, t, d, time;
cin >> s >> t >> d >> time;
s--, t--;
G[s][t].d = d;
G[s][t].time = time;
G[t][s].d = d;
G[t][s].time = time;
}
pair<ll, ll> ans = rec((1 << N) - 1, 0);
if (ans.first != INF) {
cout << ans.first << " " << ans.second << endl;
} else {
cout << "IMPOSSIBLE" << endl;
}
return 0;
}
D – 部活のスケジュール表 (Schedule)
J君・O君・I君の3人が部活に参加するかどうかのN日間のスケジュールの決め方が何通りあるか知りたい。ただし、スケジュールは以下を満たすようにしなければならない。
- 鍵を持っている人が参加しなくてはならない
- 予め定められたその日の責任者は必ず参加する
- 鍵はその日参加した人のうち誰が持ち帰っても良い
- 鍵は一つしかない
考えられるスケジュールの決め方を10007で割った余りで答えよ。
リンク：第１３回日本情報オリンピック 予選（オンライン）D – 部活のスケジュール表 (Schedule)
先程までのトラベリングセールスマンの話と少し毛色が異なりますが、これも集合を添え字にした動的計画法が利用できる問題です。つまりビットDPが使えます。
考え方
参加者の集合をビットで表現できる
1日ごとの参加者がどうなるかを考えてみると、J君・O君・I君のそれぞれについて、参加したかしなかったかの2通りで表現することができます。
つまり、J君を0bit目, O君を1bit目, I君を2bit目として参加したかどうかを0,1で表すことを考えると、参加者の集合を3bitの2進数で表すことができるのです。
例えば、1日目に J と I が参加した時を考えてみましょう。このときの参加者を2進数で表現すると \(101_2\) などと表すことができます。
動的計画法が使える問題だと見ぬく
鍵の受け渡しを考えると、前日の参加者がどうだったかによって、当日は誰が来るべきで誰が来なくても良いのかが変わってきます。
それを踏まえると、以下のように DP を定義すれば、漸化式を作れることが分かります。
dp[ 参加者の集合S ][ i ] := i 日目に参加者の集合が S だった時の、0~i 日目までのスケジュールの決め方の総和
漸化式の更新は以下のようになります。
- dp[next][i + 1] = \(\sum_{now \in \{J,O,I\}}\) dp[now][ i ] (ただし next は i+1 日目の責任者を含み、now と next で共通の参加者がいる時のみを考える)
はじめに部室の鍵を持っているのは J 君なので、初期化は以下のようにしてやれば大丈夫です。
- dp[ { J } ][ 0 ] = 1
- それ以外は全て 0
実装例
配るDPでの実装例です。実際には 10007で割った余りを出力しなければならないので忘れないようにしましょう。
#include <bits/stdc++.h>
#define LOOP(n) for (int _i = 0; _i < (n); _i++)
#define REP(i, n) for (int i = 0; i < (n); ++i)
#define RREP(i, n) for (int i = (n); i >= 0; --i)
#define FOR(i, r, n) for (int i = (r); i < (n); ++i)
#define ALL(obj) begin(obj), end(obj)
using namespace std;
using ll = long long;
const int MOD = 10007;
template <class T>
T add_mod(T &a, const T &b) {
a += b;
a %= MOD;
return a;
}
int N;
string S;
int dp[1 << 3][1100];
int resp[1100]; // resp[i]: i日目に誰が来たか？（J:1<<0, O:1<<1, I:1<<2）
int main() {
cin >> N >> S;
REP(i, N) {
if (S[i] == 'J') {
resp[i + 1] = 1 << 0;
} else if (S[i] == 'O') {
resp[i + 1] = 1 << 1;
} else if (S[i] == 'I') {
resp[i + 1] = 1 << 2;
}
}
dp[1][0] = 1; // 0日目にJ君一人が来たとする
REP(i, N) {
REP(now, 1 << 3) {
REP(next, 1 << 3) {
if (now & next) { // 2日連続で来られる人がいるか
if (next & resp[i + 1]) { // 責任者が来た時のみ加算
add_mod(dp[next][i + 1], dp[now][i]);
}
}
}
}
}
int ans = 0;
REP(now, 1 << 3) { add_mod(ans, dp[now][N]); }
cout << ans << endl;
return 0;
}
練習問題
- AtCoder Educational DP Contest O – Matching
- AtCoder Educational DP Contest U – Grouping
ビットDPを利用して計算ができます。遷移を少し工夫することで \(O(3^N)\) で計算をするタイプのものです。 - AtCoder Beginner Contest 180 E – Traveling Salesman among Aerial Cities
- AtCoder Regular Contest 016 C – ソーシャルゲーム
確率や期待値を用います。今までに入手したアイドルの集合をビットで管理するなどしましょう。