最長増加部分列(LIS)の期待値についての考察(その１：LISの導入)
こんにちは、とぼです。
最近、とても眠いです。1日に10時間くらい寝てしまう。これも春の影響でしょう か。新しい年度の始まる春ですから、もう少し気を引き締めていきたいものです。
さて、今回は動的計画法の定番問題LIS(最長増加部分列)の関する話題です。
完全に個人的な気分で書きます。メモ代わり。
でも、なるべく詳しく、分かりやすく書けるよう頑張ります。
そもそも、最長増加部分列(LIS)問題とは？
日本のオンラインジャッジサイトで最大級のAOJ（Aizu Online Judge）に典型的な最長増加部分列の問題が掲載されています。
例として、
A={ 1, 2, 0, 3, 2, 5 }
という数列を考えてみましょう。
この数列から、順番を変えずにいくつかの数字を取り出したものを「部分列」といいます。
例えば、この数列から、左から1,4,6番目の数字を取り出すと
1, 3, 5
という長さ３の部分列を得られます。この部分列は、左から右に向かって数字が大きくなっています。
つまり、1 < 3 <5 ですから (一番目の数字) < (二番目の数字) <(三番目の数字)
となっています。
このように、左から右に向かって数字が大きくなっている部分列のことを「増加部分列」といいます。
一方、この数列から、左から1,2,3,4番目の数字を取り出すと
1, 2, 0, 3
という長さ４の部分列を得られます。この部分列は、２番目の数字「2」から３番目の数字「0」にかけて数字が減ってしまっています。従って、これは増加部分列とはいえません。
また、左から1,2,4,6番目の数字を取り出すと
1, 2, 3, 5
となります。この部分列は左から右に数字が増加しているので「増加部分列」であるといえます。また、この増加部分列の長さは４であり、これより長い増加部分列は存在しないことが分かります。よって、数列Aの「最長増加部分列 = Longest Increasing Subsequence(LIS)」は「1, 2, 3, 5」であり、その長さは４であるといえます。
このように、与えられた数列の最長増加部分列の長さを求める問題を「最長増加部分列(LIS)問題」といいます。
最長増加部分列(LIS)問題を解こう！
では、最初に掲示したAOJの問題を考えてみましょう。制限時間は１秒です。「制約」を見てみると、数列の長さnは１以上100,000=10⁵以下であることが分かります。
一般的なPCやオンラインジャッジで１秒間で計算できる演算量はおおよそ10⁸回であることが知られています。ですから、長さnの数列に対しておよそn²回の演算をするとn=10⁵のとき、n²=10¹⁰（回）の演算をすることになってしまい、10⁸を優に超えてしまうので、O(n²)の（=n²に比例する演算回数の) アルゴリズムでは時間に間に合わないことが分かります。
一方で、対数log(n)を計算するとlog(10⁵)=11.51..になることが分かり、
10⁵log(10⁵)≒1.151×10⁶ < 10⁸ となるので、O(nlog(n))のアルゴリズムであれば制限時間１秒に余裕で間に合いそうだと考察できます。
つまり、LISを求めるO(nlog(n))のアルゴリズムを考えればよさそうです。
＜LISのアルゴリズムの考察＞
長さnの数列Aから数字をいくつか取り出して部分列を作る方法を2ⁿ通りあります。（n個の数字それぞれに対して、部分列の要素として「選ぶ」か「選ばない」かの２通りが考えられるからです。）
従って、全探索の計算量はO(2ⁿ)となり、全く間に合わないことが分かります(2¹⁰⁰⁰⁰⁰ ≒ 10³⁰⁰⁰⁰)。
そこで、配列を用いて計算結果をメモリに記録し、繰り返しの演算を避ける「動的計画法 = DP」を用いて解くことを考えます。
例えば、上手いDPを作って、各要素i=1～nに対して、O(n)で要素iを含むようなLISを求めたとします。しかし、このようなDPの計算量は n×O(n)=O(n²)となるので、前述の議論により結局間に合わないと考えられます。O(nlog(n))でないと間に合わないのです。そもそも計算量にlog(n)が出てくるのはどのようなときでしょうか。
結論から言うと、計算量にlog(n)が含まれるようなアルゴリズムは「二分探索」が用いられることが多いです。「二分探索」では、配列の要素を真ん中から順に条件を満たすかそうか判定していき、その結果によって次に調べる要素を上位半分の真ん中の要素にするか下位半分の真ん中の要素にするかを決定します。このようにすると毎回探索の範囲が1/2になるので、およそlog₂(n)=O(log(n))の計算量で1回の探索を終えることができます。ただし、「二分探索」は基本的に整列済みの配列（数値が大きい順や小さい順になっている）にしか適用できません。
二分探索を組み合わせることで、要素i=1~nに対してO(log(n))でLISを求める上手いDPが作れれば、全体の計算量がn×O(log(n))=O(nlog(n))となり、時間内に問題が解けそうだと分かります。
では、これを実現することで知られているDPを紹介します。なお、言語はC++です。C++の便利な関数をいくつか使っていきます。：
以下のような配列を用意します。
DP[n] : DP[i]は長さがi+1である増加部分列の最後の(一番右の)要素となりうる最小値(存在しなければINF＝便宜的に無限に大きい値)
このとき、配列DPは明らかに単調増加する（部分増加列が長くなるに連れて、増加部分列の最後の要素は確実に大きくなる）ので、数値は小さい順になるため、二分探索が適用できることが分かります。
まず、配列DPの全要素をINFで初期化します(fill関数が使える)。
次に、数列Aの1番目の要素A[0]から順に、二分探索によってA[i]以上の値をもつ配列DPの最小のインデックス（一番左の要素)を探し、そのインデックスをもつDPの要素の値をそのままA[i]に更新します(lower_bound関数で二分探索ができます)。これによって更新された(値がINFでなくなった)、
「(DPのインデックス)＋１」の長さの増加部分列は構成可能であって、
逆に数列Aの全要素に対して以上の操作を行った後も値が依然INFであるようなDPのインデックスの意味する長さの増加部分列は構成不可能であるということになりますから、「数列Aのすべての要素に対して以上の操作を行った後に値がINFでなくなったDPのインデックスの最大値(一番右にあるやつ)＋１」が数列AのLISであると分かります。
ちょっと文面だけでは分かりにくいので、具体例として最初に提示した数列
A={ 1, 2, 0, 3, 2, 5 }
に対して上記のDP配列を埋めていきましょう。
最初はDP配列の各値はINFになっています。
A[0]=1ですから、A[0]より大きい値を持つ最小のインデックスのDPはDP[0](=INF)です。よって、DP[0]=１に更新します。
次にA[1]=2ですから、A[1]より大きい値を持つ最小のインデックスのDPはDP[1](=INF)です。DP[0]の値は今１なのでA[1]=2より小さくなっていることに注意です。これより、DP[1]=2に値を更新します。
次A[2]=0ですから、A[2]より大きい値を持つ最小のインデックスのDPはDP[0](=1)です。よってDP[0]=0に更新します。
このような操作を繰り返すと、以下表のようになります。
| Aの数値 | DP[0] | DP[1] | DP[2] | DP[3] | DP[4] | DP[5] |
| (最初) | INF | INF | INF | INF | INF | INF |
| 1 | 1 | INF | INF | INF | INF | INF |
| 2 | 1 | 2 | INF | INF | INF | INF |
| 0 | 0 | 2 | INF | INF | INF | INF |
| 3 | 0 | 2 | 3 | INF | INF | INF |
| 2 | 0 | 2 | 3 | INF | INF | INF |
| 5 | 0 | 2 | 3 | 5 | INF |
INF |
この表により、最終的な値がINFでないDPの最大のインデックスは3であるので、
求めるLISは3+1=4であることが分かりました。この結果が冒頭で計算したものと確かに等しくなっていることも分かります。
以上から、このDPが正しく動作することが分かるので、あとは実装するだけです。以下のようなシンプルなコードで実装できます。
//最長増加部分列の長さの計算
int LIS(vector<int> V) {
int DP[MAX_n] = { 0 };
fill(DP, DP + V.size(), INF); //DPの各値をINFで初期化
for(int i=0;i<V.size();++i)
*lower_bound(DP, DP + V.size(), V[i]) = V[i]; //V[i]を超える最初のDPに代入
return lower_bound(DP, DP + V.size(), INF) - DP; //INFでない最大index+1
}
int LIS(vector<int> V) {
int DP[MAX_n] = { 0 };
fill(DP, DP + V.size(), INF); //DPの各値をINFで初期化
for(int i=0;i<V.size();++i)
*lower_bound(DP, DP + V.size(), V[i]) = V[i]; //V[i]を超える最初のDPに代入
return lower_bound(DP, DP + V.size(), INF) - DP; //INFでない最大index+1
}
なお、 MAX_nは数列Aの要素数として考えられる最大値を示します。先の問題ではMAX_n=100,000ですので定数constとして宣言しておけばOKです。
これでAOJの問題を解けます。
実際に今私が書いたコードを以下に示します。
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <cmath>
#include <set>
#include <map>
#include <string>
#include <numeric>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iterator>
#include <cstring>
using namespace std;
const int INF = (1<<30);
const int MAX_n = 100000;
//最長増加部分列の長さの計算
int LIS(vector<int> V) {
int DP[MAX_n] = { 0 };
fill(DP, DP + V.size(), INF);
for (int i = 0; i < V.size();++i) {
*lower_bound(DP, DP + V.size(), V[i]) = V[i];
}
return lower_bound(DP, DP + V.size(), INF) - DP;
}
int main()
{
int n;
cin >> n;
vector<int> A(n);
for (int i = 0; i < n; ++i) cin >> A[i];
cout << LIS(A) << endl;
return 0;
}
このコードでACできます。Aの要素の最大値が10⁹ということで、INF=(1<<30)、つまり2^30にしました。LIS(A)って部分、アニソン感あって良くないですか？（は）
...長くなりすぎたので、今回はLIS問題の導入ということで。あまり記事が長くなると読み込みが重くなるなどの弊害があるので。
その弊害を起こした記事👇
yadrfu.hatenablog.com ...ということで、次回、LISの期待値について考察を進めていきます！
次回が本題ではあるが、正直今回より短くなりそう（は？）。
