最長増加部分列の長さ取得アルゴリズムLISをpythonで書いてみる

最長増加部分列(Longest Increasing Subsequence )とは

数列 \{ a_i\}_{i=1}^n が与えられたとき, その部分列 \{ a_{\rho(i)}\}_{i=1}^{n'} であって
1\leqq i< j\leqq n' のとき a_{\rho(i)}< a_{\rho(j)} を満たすものを増加部分列と言います.
その中で最も長いものの長さを求めるアルゴリズムを以下に説明していきます.

1. 計算量が\mathcal{O}(n^2)アルゴリズム

まずは, 帰納的に長さを求めることを考えましょう.

\textrm{dp}(i)=\textrm{"最後が }a_i\textrm{ となるような最長増加部分列の長さ"}
とすると, 次のようなアルゴリズムが書けます.

def LIS(L):
  res=0
  dp=[1]*N
  for i in range(N):
    for j in range(i):
      if L[j]<L[i]:
        dp[i]=max(dp[i],dp[j]+1)
    res=max(res,dp[i])
  return res

i, j についての二重ループが含まれているため, このアルゴリズムの計算量は\mathcal{O}(n^2).

考え方

mn より小さい自然数とする. 仮に任意の k=1,\cdots,m に対して数列 \{ a_i\}_{i=1}^k の増加部分列のうち, 最後が a_k となるようなもののうちの最長増加部分列の長さ\textrm{dp}(k)が与えられたとしましょう.
このとき, 数列\{ a_i\}_{i=1}^{m+1}の増加部分列のうち, 最後がa_{m+1}となるような最長部分列の長さは,
\left\{
\begin{split}
    &\textrm{dp}(1)=1,\\
    &\textrm{dp}(m+1)=
      \begin{cases}
       1, & \text{任意の$k=1,\cdots,m$について$a_k \geqq a_{m+1}$のとき},\\
       \sup{\big\{\textrm{dp}(k)+1|\, a_k < a_{m+1},k=1,\cdots,m\big\}},  & \text{それ以外のとき}
      \end{cases}
\end{split}
\right.
という漸化式で与えらることを用いています. ちなみに, 漸化式の2行目ですが,
\textrm{dp}(m+1)=\max\big\{1,\sup{\big\{\textrm{dp}(k)+1|\, a_k < a_{m+1},k=1,\cdots,m\big\}\big\}}
あるいは,  A\subset \mathbb{R}_+=(0,\infty) に対して  \textrm{sup}_+

\textrm{sup}_+(A)=
  \begin{cases}
  0,& \text{$A$が空集合のとき},\\
  \sup{(A)}, & \text{それ以外のとき}
  \end{cases}
と定めたとき,
\textrm{dp}(m+1)=\sup_+{\big\{\textrm{dp}(k)|\, a_k < a_{m+1},k=1,\cdots,m\big\}}+1
と表すことができます.

2. 計算量が\mathcal{O}(n\log{n})アルゴリズム

\textrm{dp}(l)="l\textrm{ 個の増加部分列の最後の要素となる数字のうち最小の値"}
とすると, 下記のアルゴリズムで求めることが出来ます.

def LIS(L):
  from bisect import bisect
  seq=[]
  for i in L:
    pos=bisect(seq,i)
    if len(best)<=pos:
      seq.append(i)
    else:
      seq[pos]=i
  return len(seq)

i のループが含まれている一方で, 要素の追加ラベルを探すのには二分探索を用いているため, このアルゴリズムの計算量は\mathcal{O}(n\log{n}). ただし seq そのものは最長増加部分列になっているわけではないので注意しましょう.

考え方

 \textrm{dp} の長さを求めます.
 l を固定して, \textrm{dp}(l) の値を数列 \{a_j\}_{j=1}^n の第 i 項までを用いて表すことを考えると, これは i について単調減少になるはずなので, その場合に更新することを考えています.
更新する際に考えなくてはならないのは a_k \textrm{dp} の第何項に該当するかですが,  \textrm{dp} は常に単調減少な配列になっていることから二分探索 bisect を用いることができて計算量を減らしています.

次回はこれを用いた応用として競技プログラミングの問題を考えていきます.