最長増加部分列の長さ取得アルゴリズム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{seq}(l)="\textrm{長さ }l\textrm{ の増加部分列の最後の要素となりうる数字のうち最小の値"}
とすると, 下記のアルゴリズムで求めることが出来ます.

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

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

考え方

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

3. 広義単調増大部分列の場合

追加された要素の挿入位置を探す際に, 二分探索して同じ数値が得られる場合には考えられる右端の index を返すようにします.
このためには, `bisect_left` の代わりに, `bisect` もしくは `bisect_right` を用います.
bisect --- 配列二分法アルゴリズム — Python 3.9.4 ドキュメント

def improper_LIS(L):
    from bisect import bisect
    seq = []
    for ai in L:
        pos = bisect(seq, ai)
        if len(seq) <= pos:
            seq.append(ai)
        else:
            seq[pos] = ai
    return len(seq)
次回はこれを用いた応用として競技プログラミングの問題を考えていきます.