「競プロ典型 90問」Smallest Subsequence (最小部分列問題)
最小部分列問題
「 競プロ典型 90 問」の 006 - Smallest Subsequence（★5） (最少部分列問題) という問題を解いてみたのですが、最初は解説をみてもさっぱり分からず打ちひしがれていました･･･。
が、けんちょんの競プロ精進記録 を見るに、どうもこの問題を解く途中で出てくる nex という配列が「極めて汎用性が高いので、実にさまざまな問題で活用できます！！！」ということらしく、ちゃんと理解しといた方が良さそうだ･･･ということで気を取り直して取り組んでみたところなんとか理解できました。
せっかくなので忘れないうちに解説記事を作って記憶を定着させたいと思います。なお後半の実装パートは、Haskell で実装します。
けんちょんさんの解説記事にあるとおり、この問題 (を全探索で解く場合) の解法のキーになるのは事前に「任意の文字が
最小部分列とは何か?
- 長さ
N S - 長さが
K S
「部分列」は要するに、
kittyonyourlap
の 5 文字の(最小ではない) 部分列は順序さえ変えなければいいので
kittyonyourlap -> kyoyo
* ** **
も5文字の部分列。こうやって作られる部分列のうち、辞書式順で最小のものが最小部分列。kittyonyourlap の場合は
kittyonyourlap -> inlap
* * ***
この inlap が最小部分列になります。
ざっとみたときに、たしかに5文字の部分列を作ろうとすると
k*****
i***** <- k 始まりよりも i 始まりのほうが辞書式順で小さい
t***** <- i 始まりの方が小さい
...
r*** <- r が先頭になると5文字作れない
l** <- 以降 5文字が作れない
a*
p
となるので、確かに i を先頭にした部分列が最小部分列になるのが分かります。
i を先頭にした部分列は
kittyonyourlap -> ittyo
*****
など他にも幾つか考えられるのですが、最小部分列となると ittyo などよりも inlap の方が辞書式順では小さい、ということになります。
最小部分列の求め方の基本的な考え方
最小部分列を求めていく基本的な考え方としては、途中でみた
k*****
i***** <- k 始まりよりも i 始まりのほうが辞書式順で小さい
t***** <- i 始まりの方が小さい
...
r*** <- r が先頭になると5文字作れない
l** <- 以降 5文字が作れない
a*
p
この考え方がもとになります。
kittyonyourlap
****
後ろの kittyonyour のうち一番小さな文字である
1文字目は i で確定したので、2文字目を考えます。5文字中1文字が確定したので、次は 4文字の部分列を考えればよいです。このとき、1文字目は
k
i <- 1文字目
--- 2文字目はここから先しか考えなくて良い --
t***
t***
y***
o***
n*** <- これ
y***
o***
u***
r***
--- ここ以降は検討しなくてよい --
l** <- 4文字にならない
a*
p
となります。
結果
もう1文字やってみましょう。3文字目です。
y**
o**
u**
r**
l** <- これ
a* <-- 対象外
p <-- 対象外
(i, 2) -> (n, 7) -> (l, 12) -> (a, 13) -> (p, 14)
という流れで inlap という最小部分列を求めることができます。
整理すると kittyonyourlap という
- 5文字を構成できるのは
14 - (5 - 1) = 10 S[1] \space .. \space S[10] S[2] = i - 2文字目は
S[2] = i S[3] \space .. \space S[11] S[7] = n - 3文字目は
S[7] = n S[8] \space .. \space S[12] S[12] = l - 4文字目は
S[12] =l S[13] \space .. \space S[13] S[13] = a - 5文字目は ･･･
S[14] = p
となります。
こうしてみてみると「
a 〜 z を全探索で探す方法
「a から順にその範囲をみていって、a が含まれてなかったら次は b 、その次は c ･･･とみていって、見つかった文字 (とその位置) を返すという愚直な方法でも十分に高速です。
例えば 5文字を構成できる先頭文字を探す場合
1 2 3 4 5 6 7 8 9 10
----------------------------
k i t t y o n y o u
の範囲を a から順に当てていくと i がインデックス 2 で見つかります。
次は、この i の位置 a 〜 z を当てていきます。
3 4 5 6 7 8 9 10 11
-------------------------------
t t y o n y o u r
a b c d e ... となかなかヒットしないですが、そのうち n が位置
こうして考えると 「文字列
「i c
この「文字列
例えば kittyonyourlap に対する nex と配列が既に構築されていたとして
print $ nex ! 1 ! 'i' -- 2
print $ nex ! 3 ! 'n' -- 7
print $ nex ! 8 ! 'l' -- 12
print $ nex ! 13 ! 'a' -- 13
print $ nex ! 14 ! 'p' -- 14
-- みつからない場合
print $ nex ! 1 ! 'z' -- -1
という具合に使うことができる配列。これをあらかじめ構築します。
で、今回なかなか理解に至らなかったのが実はこの配列の構築方法でして「文字列
構築方法よりも先に、まずできあがりの配列がどんな形をしているのかを考えてみます。
a b c d e f g h i j k l m n o p q r s t u v w x y z
-----------------------------------------------------------------------------------
1 k 13 2 1 12 7 6 14 11 3 5
2 i 13 2 12 7 6 14 11 3 5
3 t 13 12 7 6 14 11 3 5
4 t 13 12 7 6 14 11 4 5
...
こういう配列になっています。
例えば inlap の i を決めるのに a 列からみていくのが a .. z を当てていくことに相当します。すると最初に a = 13 が見つかるのですが、i = 2 が見つかる。
inlap の n を決めるのには a = 13 と l = 12 が先にみつかりますがいずれも範囲外で、n = 7 に当たります。
･･･というのが、欲しい nex 配列です。
これをどう構築するかですが、「文字列の
a b c d e f g h i j k l m n o p q r s t u v w x y z
-----------------------------------------------------------------------------------
14 p 14
13 a 13 14
12 l 13 12 14
11 r 13 12 14 11
10 u 13 12 14 11
9 o 13 12 9 14 11
8 y 13 12 9 14 11 8
7 n 13 12 7 9 14 11 8
6 o 13 12 7 6 14 11 8
5 y 13 12 7 6 14 11 5
4 t 13 12 7 6 14 11 4 5
3 t 13 12 7 6 14 11 3 5
2 i 13 2 12 7 6 14 11 3 5
1 k 13 2 1 12 7 6 14 11 3 5
ふむ。これで構築方法がぱっと思いつくか。難しいところですが、最初の数行に絞ってみてみます。
a b c d e f g h i j k l m n o p q r s t u v w x y z
-----------------------------------------------------------------------------------
14 p 14
行の p しか出現せず、その位置は (行のインデックス, 列の文字) としたとき、最初に (14, p) = 14 が求まります。
a b c d e f g h i j k l m n o p q r s t u v w x y z
-----------------------------------------------------------------------------------
14 p 14
13 a 13 14
その次は (13, a) = 13 が apであり、そこには p も出現します。p の出現位置は、当然先に同じ (13, p) = 14 も埋めます。これは一つ前の行からそのまま値を次の行に下ろしてやる (コピーしてやる) ので良いでしょう。
a b c d e f g h i j k l m n o p q r s t u v w x y z
-----------------------------------------------------------------------------------
14 p 14
13 a 13 14
12 l 13 12 14
次、(12, l) = 12 が決まります。一方で、a や p はやはり出現位置が同じなので上から下にコピーする。
こうして 1行遷移するごとに、
･･･二次元配列を構成するのに、一つ前の行の計算結果を、次の行の計算結果に使う。この遷移を繰り返す。なんだかどこかで見たようなパターンですね。そうです、二次元DPのDP表を構築するのと同じような計算です。
実際に、DP とは違って緩和 (複数のセルの値が特定のセルに合流するケース) があるわけではないので、DP ではなくただの累積計算なのですが二次元表を構築していくフレームワークとして二次元DPのそれと同じように考えると (二次元DP に慣れていれば) すんなり理解できます!
# res[i][c] := i 文字目以降で最初に文字 c が登場する index
# 存在しないときは N
def calc_next(S):
# 文字列 S の長さ
N = len(S)
# 答え
res = [[N] * 26 for _ in range(N + 1)]
# 後ろから見る
for i in range(N - 1, -1, -1):
# i + 1 文字目以降の結果を一旦 i 文字にコピー
for j in range(26):
res[i][j] = res[i + 1][j]
# i 文字目の情報を反映させる
res[i][ord(S[i]) - ord('a')] = i
# 答えを返す
return res
けんちょんさんの解説記事から引用した Python での nex 配列の構築の実装ですが「# i + 1 文字目以降の結果を一旦 i 文字にコピー」が、先の表でいうところの前の行の結果を次の行にコピーすることに相当していたのでした。
nex 配列構築を Haskell で実装する
というわけでここから実装パートです。Haskell の話題に移っていきます。
二次元DP と同じ考え方だとわかれば、普段二次元DPを実装しているとき同様の方法に落とし込めるので簡単です。
Haskell では二次元配列の構成方法は例えば Array の Array を作る方法などもありますが、Array のインデックスにはタプルが使えて、それで二次元表を表現することもできます。
その方法で nex 配列を構築するのが以下の実装です。
import Data.Array.IArray (Array, array, (!))
import qualified Data.IntMap.Strict as IM
buildNex :: Int -> String -> Array (Int, Char) (Maybe Int)
buildNex n str = nex
where
-- 入力 s を「文字の位置 => 文字」で検索できるように Map にしておく
s :: IM.IntMap Char
s = IM.fromList (zip [1 ..] str)
nex :: Array (Int, Char) (Maybe Int)
nex = array ((1, 'a'), (succ n, 'z')) $ row0 ++ rows
-- 番兵になる最初の行。すべて Nothing
row0 :: [((Int, Char), Maybe Int)]
row0 = [((succ n, c), Nothing) | c <- ['a' .. 'z']]
-- n 行から n - 1 行･･･と降順に遷移させていく
rows :: [((Int, Char), Maybe Int)]
rows = [((i, c), f i c) | i <- [n, n - 1 .. 1], c <- ['a' .. 'z']]
where
f :: Int -> Char -> Maybe Int
f i c
| c == (s IM.! i) = Just i -- c が s に含まれているならその位置を使う
| otherwise = nex ! (i + 1, c) -- 一つ上の行の値をコピー (遅延評価で構築中の nex 配列を参照する)
main :: IO ()
main = do
let nex = buildNex 14 "kittyonyourlap"
print $ nex ! (1, 'a') -- Just 13
print $ nex ! (1, 'b') -- Nothing
print $ nex ! (1, 'i') -- Just 2
print $ nex ! (3, 'n') -- Just 7
Haskell なら遅延評価を使うことで、 nex 配列を構築する最中に nex 配列にアクセスすることができるので、それを利用して Array を構築します。DP 表を作るときにも、筆者の場合、同様の実装でまずやってみることが多いです。
これで nex 配列が期待通りに構築できました。
nex 配列を使って最小部分列を復元する
nex 配列さえ構築できればあとは簡単です。
- 前半でみたとおり、
S S[1] \space .. \space S[N - (K - 1)] a .. zの全探索。(以下、それをやっているのがnextChar関数です。) - そこでみつかった索引を
i S[i + 1] \space .. \space S[N - (K - 2)]
これを繰り返します。
import Data.Array.IArray (Array, array, (!))
import qualified Data.IntMap.Strict as IM
import Data.List (mapAccumL)
import Data.Maybe (mapMaybe)
buildNex :: Int -> String -> Array (Int, Char) (Maybe Int)
buildNex n str = nex
where
s :: IM.IntMap Char
s = IM.fromList (zip [1 ..] str)
nex :: Array (Int, Char) (Maybe Int)
nex = array ((1, 'a'), (succ n, 'z')) $ row0 ++ rows
row0 :: [((Int, Char), Maybe Int)]
row0 = [((succ n, c), Nothing) | c <- ['a' .. 'z']]
rows :: [((Int, Char), Maybe Int)]
rows = [((i, c), f i c) | i <- [n, n - 1 .. 1], c <- ['a' .. 'z']]
where
f :: Int -> Char -> Maybe Int
f i c
| c == (s IM.! i) = Just i
| otherwise = nex ! (i + 1, c)
nextChar :: Array (Int, Char) (Maybe Int) -> (Int, Int) -> (Int, Char)
nextChar nex (s, e) = head $ mapMaybe f ['a' .. 'z']
where
f :: Char -> Maybe (Int, Char)
f c = do
i <- nex ! (s, c)
if i <= e then return (i, c) else Nothing
solve :: Int -> Int -> String -> String
solve n k s = snd $ mapAccumL f 0 [k - 1, k - 2 .. 0]
where
nex = buildNex n s
f pos l = nextChar nex (succ pos, n - l)
main :: IO ()
main = do
[n, k] <- map read . words <$> getLine
s <- getLine
putStrLn $ solve n k s
これで無事 AC しました。
･･･が、メモリを 200MB ぐらい使った上に実行時間が 604 ms 。C++ で AC している人たちの解答をみると 40 ms 弱〜 200ms ぐらい。Python 3 でも 100 〜 400ms ぐらいですので、この Haskell の実装は遅いですね。
遅延評価のために Array を使ってサンクを積みまくる実装なので、遅いのも仕方がない。プロトタイプ実装みたいなもんだと思ってください。
STUArray で高速にする
というわけで、改善。遅延評価は捨てて STUArray を使います。実装の方針はほぼ一緒です。
import Control.Monad (forM_)
import Control.Monad.ST (ST)
import Data.Array.IArray ((!))
import Data.Array.MArray (
MArray (newArray),
readArray,
writeArray,
)
import Data.Array.ST (STUArray, runSTUArray)
import Data.Array.Unboxed (UArray)
import qualified Data.ByteString.Char8 as BS
import qualified Data.IntMap.Strict as IM
import Data.List (mapAccumL)
import Data.Maybe (mapMaybe)
buildNex :: Int -> String -> UArray (Int, Char) Int
buildNex n str = runSTUArray $ do
nex <- newArray ((1, 'a'), (n + 1, 'z')) (-1)
forM_ [(i, c) | i <- [n, n - 1 .. 1], c <- ['a' .. 'z']] $
\x -> f nex x >>= writeArray nex x
return nex
where
s :: IM.IntMap Char
s = IM.fromList (zip [1 ..] str)
f :: STUArray s (Int, Char) Int -> (Int, Char) -> ST s Int
f nex (i, c)
| c == (s IM.! i) = return i
| otherwise = readArray nex (i + 1, c)
nextChar :: UArray (Int, Char) Int -> (Int, Int) -> (Int, Char)
nextChar nex (s, e) = head $ mapMaybe f ['a' .. 'z']
where
f :: Char -> Maybe (Int, Char)
f c = do
let i = nex ! (s, c)
if i > (-1) && i <= e then return (i, c) else Nothing
solve :: Int -> Int -> String -> String
solve n k s = snd $ mapAccumL f 0 [k - 1, k - 2 .. 0]
where
nex = buildNex n s
f pos l = nextChar nex (succ pos, n - l)
main :: IO ()
main = do
[n, k] <- map read . words <$> getLine
s <- BS.unpack <$> BS.getLine
putStrLn $ solve n k s
これで 74ms 。だいぶ速くなりました。悪くないでしょう。
イミュータブルに nex 配列を構築するより良い方法
高速化するのに STUArray を使いましたが joetheshootingst さんの解答 を参考に、イミュータブルな Array でももっとエレガントな実装にしてみたのが以下です。
配列の各行の遷移 ･･･ 状態の遷移を scanr で実装することでサンクを積み過ぎることなく nex 配列をイミュータブルに構築することができます。nex 配列の各行は、一つ前の行からの遷移を表すわけですがこの問題の場合その遷移元は必ず一つ前の行のみになることが分かっているので scanr がはまるわけですね。
{-# LANGUAGE TypeApplications #-}
import Data.Array.IArray (Array, accumArray, listArray, (!), (//))
import Data.Array.Unboxed (UArray)
import qualified Data.ByteString.Char8 as BS
import Data.List (mapAccumL)
import Data.Maybe (mapMaybe)
type Nex = Array Int (UArray Char Int)
buildNex :: Int -> String -> Nex
buildNex n s = listArray @Array (1, n) rows
where
row0 = accumArray @UArray const (-1) ('a', 'z') [] -- 初期値の埋まった配列を accumArray で構築して･･･
rows = scanr (\(c, i) pos -> pos // [(c, i)]) row0 $ zip s [1 ..] -- それをベースに scanr で遷移させながら、段階的に配列を更新
nextChar :: Nex -> (Int, Int) -> (Int, Char)
nextChar nex (s, e) = head $ mapMaybe f ['a' .. 'z']
where
f :: Char -> Maybe (Int, Char)
f c = do
let i = nex ! s ! c
if i > (-1) && i <= e then return (i, c) else Nothing
solve :: Int -> Int -> String -> String
solve n k s = snd $ mapAccumL f 0 [k - 1, k - 2 .. 0]
where
nex = buildNex n s
f pos l = nextChar nex (succ pos, n - l)
main :: IO ()
main = do
[n, k] <- map read . words <$> getLine
s <- BS.unpack <$> BS.getLine
putStrLn $ solve n k s
この実装だと 123ms でした。STUArray を使った実装には及びませんが、最初の実装よりはかなり速い。とてもバランスの良い実装です。
全探索ではなくヒープを使う
ここまでは全探索で解く方法を解説してきましたが、この問題はヒープやセグメント木を使って解くこともできます。
以下がヒープを使った実装です。Haskell の Data.Heap を使います。(Data.Heap の使い方については Haskell で、優先度付きキューを使ったダイクストラ法 の中でも解説していますので、興味のある方はご一読ください)
{-# LANGUAGE TypeApplications #-}
import Data.Array.Unboxed (UArray, listArray, (!))
import qualified Data.ByteString.Char8 as BS
import qualified Data.Heap as Heap
import Data.List (mapAccumL)
import Data.Maybe (fromJust)
type Queue = Heap.Heap (Char, Int)
nextChar :: Queue -> Int -> ((Char, Int), Queue)
nextChar h start
| i >= start = ((c, i), h')
| otherwise = nextChar h' start
where
((c, i), h') = (fromJust . Heap.uncons) h
solve :: Int -> Int -> String -> String
solve n k s = snd $ mapAccumL f (0, q0) [k - 1, k - 2 .. 0]
where
sa = listArray @UArray (1, n) s
q0 = Heap.fromList $ zip s [1 .. n - k]
f :: (Int, Queue) -> Int -> ((Int, Queue), Char)
f (i, q) l = ((i', q'), c)
where
qq = Heap.insert (sa ! (n - l), n - l) q
((c, i'), q') = nextChar qq (succ i)
main :: IO ()
main = do
[n, k] <- map read . words <$> getLine
s <- BS.unpack <$> BS.getLine
putStrLn $ solve n k s
この実装で 93ms でした。
こちらの実装の方が、ロジックとしてはシンプルで認知負荷が低いかもしれないですね。
長くなってきたので解説はざっくりですが
kittyonyourlap
を例に考えたとき
- ステップ1 ･･･
*****の5文字の先頭文字を見つけたい。逆に言えばそこが先頭になっても5文字作れないrlapは探索の対象外- ヒープに
S[1] S[n - (k - 1)] - ヒープから値を取り出すと辞書式順で最少値である
(i, 2)が取り出される
- ヒープに
- ステップ2 ･･･
****の4文字の先頭文字を見つけたい。ここで先にヒープに入れなかったrlapのうちrは探索の対象になる- ヒープに
(r, 11)を insert する - ヒープから値を取り出すと
(k, 1)が取り出されるが、1 は(i, 2)の 2 よりも前なのでスキップしてまたヒープから取り出す- これを繰り返してインデックスが 3 以上の値が出てくるまでヒープから取り出す ･･･
(n, 7)が手に入る
- これを繰り返してインデックスが 3 以上の値が出てくるまでヒープから取り出す ･･･
- ヒープに
- ステップ3 ･･･
***の3文字の先頭文字を見つけたい。ヒープにまだ入っていないlapのうちlが探索対象になる- ヒープに
(l, 12)を insert する -
(n, 7)の次以降なのでインデックスが 8 以上の値が出てくるまでヒープから取り出す ･･･(l, 12)が手にはいる
- ヒープに
これを繰り返すと (i, 2) -> (n, 7) -> (l, 12) -> (a, 13) -> (p, 14) の順に pop される
-
K K - 1 - ヒープから値を取り出すと辞書式順にソートされているし、対象外になる文字はそもそもヒープに入れてないので、取り出された値が候補になる
- ただし、その取り出した値が先に取り出した文字よりも左にあるなら、スキップ
このループになる･･･ということをやっています。
そもそもなぜヒープでうまく実装できるかというと、そもそも最小部分列を求めるアルゴリズムが
- 最小部分列を先頭から文字を決めていくにあたっては、
i N - (K - j) - 一方でステップごとにこの
i j
という貪欲法に基づいているからで、この「ステップ毎にその時点での最善手を探す」という操作に、その性質上ヒープがうまく適合するわけです。この問題に限らず貪欲法ではヒープを使う場面がよくあります。
Wikipedia によれば貪欲法は
ステップからなる意志決定問題において、ある特定のルールに基づいて「その場での最善の選択」を繰り返していく方法論 N
ということですが、つまり全体を見渡しての最善の選択ではなく、ステップを繰り返す最中、あくまでその時点での最善手を考えるからでしょうね。全体を見渡しての最善手だと、予めソートされているとか、DP のような形であらゆる分岐を考慮した上で計算するということが必要になるわけですが、貪欲法はあくまで「その時点での最善手」でありDP のように分岐もないので、ただひとつのヒープのように動的に構成できるデータ構造で、この「その時点での最善手が定まる」ということがよくあるのだと思います。
本問は貪欲法の典型問題ということですが、貪欲法を実装するという観点にフォーカスして考えるとヒープを使った実装のほうがより貪欲法が適用できる問題の構造を活かした解法ではないか、という気がします。
というわけで長くなりましたが、最小部分列問題の解説でした。
2023.7.24 追記
記事執筆から半年がすぎまして、私の Haskell 力も多少なりとも向上しもう少し簡潔に記述できたのでそのコードを載せておきます。
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeApplications #-}
import Data.Array.Unboxed (UArray, accumArray, (!))
import qualified Data.ByteString.Char8 as BS
import Data.Char (isSpace)
import Data.List (mapAccumL, scanl', unfoldr)
import Data.List.Extra (minimumOn)
import qualified Data.Map.Strict as Map
-- i 文字目以降で c が登場する最初の位置はどこか? という配列 nex を構築
buildNex :: Int -> String -> UArray (Int, Char) Int
buildNex n s = accumArray max (-1) ((0, 'a'), (n + 1, 'z')) $ do
(i, state) <- scanl' (\(i, state) c -> (i - 1, Map.insert c i state)) (n, Map.empty) $ reverse s
[((i, c), j) | (c, j) <- Map.toList state]
-- nex を使って指定された LR 区間内で辞書順最小の文字とその位置を返す
findNext :: UArray (Int, Char) Int -> Int -> Int -> (Char, Int)
findNext nex l r = do
let xs = filter ((/= -1) . snd) $ map (\c -> (c, nex ! (l, c))) ['a' .. 'z']
minimumOn fst $ filter ((<= r) . snd) xs
main :: IO ()
main = do
[n, k] <- getInts
s <- getLine
let nex = buildNex n s
-- LR区間を決めてその範囲を nex で検索、を繰り返す
result = mapAccumL f (0, n - k + 1) [1 .. k]
where
f (l, r) _ = do
-- 次の文字を見つけてLR区間を更新する
let (c, i) = findNext nex l r
((i, r + 1), c)
putStrLn (snd result)
{-- Library --}
getInts :: IO [Int]
getInts = unfoldr (BS.readInt . BS.dropWhile isSpace) <$> BS.getLine
Discussion