Python|最長増加部分列の長さと、部分列の１つを再構成する方法

tadanori

最長増加部分列（LIS）とは、数列Aの部分増加列の中で最も長いものです。この記事では、これを探索するアルゴリズムを紹介します。

Contents

最長増加部分列（LIS）とは
最長増加部分列の長さを求める方法
最長増加部分列の１つを再構成する方法
まとめ

最長増加部分列（LIS）とは

最長増加部分列（Longest Increasing Subsequence, LIS）とは、数列$A_n$の部分列として作ることができる増加部分列の中で最も長い部分列のことです。

例えば、数列A=[3,1,4,2,5]の場合、長さ２以上の増加部分列として、[3,4],[3,5],[3,4,5],[1,4],[1,4,5],[1,2],[1,2,4],[4,5]などを作ることが可能です。この中で最も長いのは[3,4,5],[1,4,5],[1,2,5]で、最長増加部分列の大きさは３となります。

この記事では、長さを求めるアルゴリズムを紹介し、また、最長部分列を１つ作成する方法を紹介します。

最長増加部分列の長さを求める方法

以下が、最長増加部分列の長さを求めるプログラムになります。

import random
import bisect

def lis(a, inf = 1e18) :
    n = len(a)
    dp = [inf] * n
    for i in range(n) :
        pos = bisect.bisect_left(dp, a[i])
        dp[pos] = a[i]

    return bisect.bisect_left(dp, inf)


a = [random.randint(1, 10) for _ in range(5)]
print(a)
print(lis(a))

# 実行結果
# [5, 9, 5, 1, 9]
# 2

このプログラムは動的計画法と二分探索法を組み合わせたものになります。dp[idx]という変数には、「i個目までの値を見た時に、idx番目に小さい変数の値」が格納されています。

A=[5,9,5,1,9]で考えてみます

最初dpは全て最大値で初期化されています（[inf, inf, inf, inf, inf]）。

まずA[0] = 5で考えると、これはdpの中で一番小さい値です。ですから、dp[0]の値を更新します。

[5, inf, inf, inf, inf]

次に、a[1] = 9だとこれはdpの中で５の次に小さい値になります。なので、5の隣を9に書き換えます。

[5, 9, inf, inf, inf]

次にa[2] = 5ですが、これは1番目の値と同じです。なので５を書き換えます（同じ値なので更新しても変化しません）

[5, 9, inf, inf, inf]

次にA[3] = 1で、これは5よりも小さくなります。ここでは、5が1に書き換えられます。

[1, 9, inf, inf, inf]

最後にA[4] = 9ですが、9はすでにありますので、9の場所が更新されます（変化しない）。

[1, 9, inf, inf, inf]

最終的に、infの手前に２の値が入っているので、最長部分増加列の長さは２となります。

なお、「それぞれの値が何番目に小さいか」を判定する処理を二分探索で行います。

先頭からN回二分探索を行うので、計算量はO（NlogN）となります。

なぜこれでうまくいくのか？

dp[idx]に記録されている値は、先頭からi番目までに現れた数値でidx番目に小さい値です。つまり、現在の位置iより手前に必ずdpに記録されている値がidx番目に存在することになります。なので、idx個の増加部分列を必ず作ることができます。

例えばdp=[1,4,5,7,inf,inf]の場合、a[i]=6であれば、i番目より手前に1,4,5の３つは確実に存在していることになります。なのでi番目までの最長増加部分列は4になります。

最長増加部分列の１つを再構成する方法

最長増加部分列の長さが求まりました。次は、最長増加部分列の１つを求めてみたいと思います。最初に説明したように、数列A=[3,1,4,2,5]の場合、最長増加部分列の長さは３で、[3,4,5],[1,4,5],[1,2,5]の３パターンが存在します。ここでは、この中の１つを再構成することを考えます。

実際に作成したプログラムは以下になります。lis2()関数は、最長増加部分列の長さと、最長増加部分列の要素となる要素番号と、値を返す関数になります。

import random
import bisect

def lis2(a, inf = 1e18) :
    n = len(a)
    dp = [inf] * n
    l = []
    for i in range(n) :
        pos = bisect.bisect_left(dp, a[i])
        dp[pos] = a[i]
        l.append(pos+1)

    ans = []
    pos = []
    cur = bisect.bisect_left(dp, inf)
    for i in range(n-1, -1, -1) :
        if l[i] == cur :
            ans.append(a[i])
            pos.append(i+1)
            cur -= 1

    return bisect.bisect_left(dp, inf), pos[::-1], ans[::-1]

a = [random.randint(1, 10) for _ in range(10)]
print(a)
n, p, v =  lis2(a)
print(n)
print(*p)
print(*v)

# 実行結果
# [7, 9, 9, 7, 6, 10, 10, 8, 2, 1]
# 3
# 1 3 7
# 7 9 10

処理としては、動的計画法でdpを更新すると同時に、i番目の数値で作れる最大の長さをlに記憶しています。

lには、i番目の値で作れる最大の長さが記録されているので、それを逆順（末尾）に調べてl[i]=最長増加部分列の長さになるものを調べ、そこから長さを-1しながら一致するものを順番に探していきます。例えば、長さが3の場合、3, 2, 1という数字をlの末尾から検索します。見つかったものが最長増加部分列のパターンの1つになります。

このように再構成することで、パターンの１つを見つけることができます。