Monotone Minima と SMAWK

以下のような問題があるとする。

• 例題1
  • $N \times M$ の2次元配列 $A$ がある
  • 各行毎に、その行の最小値を求めたい

A 
i\j 0 1 2 3
 0  3 1 4 1  → 1
 1  5 9 2 6  → 2
 2  5 3 5 8  → 3
 3  9 7 9 3  → 3

これは、通常なら、各行毎にチェックして $O(NM)$ かける他にない。

ところが、$A$ がある“都合の良い条件”を満たしていると、$O((N+M) \log{N})$ や $O(N+M)$ でできる。

• Totally Monotone Matrix Searching (SMAWK algorithm) - 週刊 spaghetti_source - TopCoder部
• SMAWK Algorithm | Library

「$i$ の増加に伴って、最小値を取るindex $j$ が単調増加」の時、$A$ は「Monotone」という。

Monotone
i\j 0  1  2  3  4
 0  9 [3] 4  8  5    $i$ の増加につれ、最小値をとる $j$ は右に
 1  5 [1] 9  2  6
 2  5  7  6 [4] 8
 3  8  7  9  3 [2]
 4  7  9  8  8 [6]

また、$A$ から行と列を抜き出した任意の部分行列 $A'$ がMonotoneであれば、$A$ は「Totally Monotone」という。

Totally Monotone
i\j  0   1   2   3   4
 0  [7] 11  20  24  50
 1  23 [13] 17  16  36
 2  47  22  21 [16] 31
 3  75  42  37 [25] 34

↓たとえば、i=[0,2,3], j=[1,2,4] を抜き出してみても...

i\j  1   2   4    部分行列でも
 0 [11] 20  50    Monotone性が保たれている
 2  22 [21] 31
 3  42  37 [34]

元の例題1に立ち返ると、$A$ が、

• Monotone であれば、Monotone Minima により、$O(N+M \log{N})$
• Totally Monotone であれば、SMAWK により、$O(N+M)$

で各行の最小値を求める方法がある。

実際に競技プログラミングの問題として Monotone 性を利用するのが出てくる場合、 $A$ は何かしら規則的な方法で生成されるのが大半だろう。

その場合「Monge性」を持つことを使うことが多い（らしい）。
「Monge性を満たすならTotally Monotone性も満たす」という性質がある。

• [Monge まとめ①] Monge 性とは？ - HackMD

• monotone minima による高速化について | フューチャー技術ブログ

$i$ について半分ずつ範囲を狭めていく。

i\j  0 1 2 3 4 5 6 7 8 9    N=7  M=10 として
 0               x x x x
 1               x x x x    まず、Nの真ん中の i=3 について O(M) かけて最小値を探す
 2               x x x x
 3   # # # # #[#]# # # #    j=5 のときに最小値を取るとすると、
 4   x x x x x              Monotone性より、左下と右上は可能性を除外できる。
 5   x x x x x
 6   x x x x x              結果、左上と右上の2つの矩形が残る。

※ #:値を調べた箇所  [#]:調べた中での最小値  x:最小値の可能性がなくなった箇所


i\j  0 1 2 3 4 5 6 7 8 9    分割された2つの矩形で再帰的に同じことをする。
 0       x x x x - - - -    行 [0,1,2] の真ん中 i=1 と、
 1   #[#]# # # # - - - -    行 [4,5,6] の真ん中 i=5 について探索する
 2   x           - - - -    この探索幅は、1箇所のみ重複するが、あわせておよそ O(M) である。
 3   - - - - - - - - - -    
 4   - - - - -       x x    最小値を見つけたら、
 5   - - - - - # #[#]# #    同様に、それぞれの左下と右上は可能性を除外でき、4つの矩形が残る。
 6   - - - - - x x          
 
 ※ -:前回までで調べた or 可能性がなくなり、調べる必要が無くなった箇所


i\j  0 1 2 3 4 5 6 7 8 9    再帰的に同じことをする。
 0  [#]# - - - - - - - -    探索幅は、あわせて O(M) である。
 1   - - - - - - - - - -    
 2   - # #[#]# # - - - -    
 3   - - - - - - - - - -    log N 回繰り返すと、全ての i について最小値の位置が判明する。
 4   - - - - -[#]# # - -    
 5   - - - - - - - - - -
 6   - - - - - - - # #[#]

探索幅は $O(M)$ と言っているが、厳密には重複が発生している。その回数は、全体を通して $O(N)$ となる。よって、全体で $O(N+M \log{N})$ で求められる。

• SMAWK Algorithm | Library
• わかりやすい具体例が載っている論文（英語, pdf）
  • The SMAWK Algorithm

定数倍の関係で、必ずしも Monotone Minima より早いとは限らない？

• GitHub - ngtkana/bench_minplus_convolutions: Benchmarks of algorithms for convex min-plus convolutions

Monotone Minima の方がわかりやすく、実装も容易で、計算量の差もlog1つだけなので、 競技プログラミングでは無理してこちらを使う必要性は薄いかも？

アルゴリズムの日本語の文章化は以下のブログが詳しい。

• Totally Monotone Matrix Searching (SMAWK algorithm) - 週刊 spaghetti_source - TopCoder部

Python3 による実装

def smawk(n: int, m: int, select: Callable[[int, int, int], bool]) -> List[int]:
    """
    n×m 配列 A の、各行における最小値(※)を与える列index [j1,j2,...,jn] を得る。

    制約:
      A はTotally Monotone
      A は陽に与える必要は無く、比較関数 select によって評価する
        select(i,j,k) := A[i,j] > A[i,k] なら True を返す関数

    (※)
    最小値以外（最大値など）の指標にも適用できる。
    Aはその指標にとっての Totally Monotone であればよく、select(i,j,k) はA[i,k]の方が"よい"場合にTrueを返す関数とする。

    参考:
      https://github.com/noshi91/Library/blob/master/algorithm/smawk.cpp
    """
    columns_history = [list(range(m))]
    row_bit_length = n.bit_length()
    for depth in range(row_bit_length):
        next_row_size = n >> depth
        i_step = 1 << depth
        i = i_step - 1
        columns = []
        for j in columns_history[depth]:
            while columns and select(i - i_step, columns[-1], j):
                columns.pop()
                i -= i_step
            if len(columns) < next_row_size:
                columns.append(j)
                i += i_step
        columns_history.append(columns)

    assert len(columns_history[row_bit_length]) == 1

    result = [0] * n
    result[(1 << (row_bit_length - 1)) - 1] = columns_history[row_bit_length][0]

    for depth in range(row_bit_length - 1, 0, -1):
        i_step = 1 << depth
        h_step = i_step >> 1
        columns = columns_history[depth]
        ji = 0
        for i in range(h_step - 1, n, i_step):
            res = columns[ji]
            until = columns[-1] if i + h_step >= n else result[i + h_step]
            while columns[ji] < until:
                ji += 1
                if select(i, res, columns[ji]):
                    res = columns[ji]
            result[i] = res

    return result

Monotone Minima や SMAWK を使うと、$C_i=\min_j(A_j+B_{i-j})$ みたいな $C$ も、$A,B$ の性質次第で高速に求められる。

• monotone や monge 周りの話 - PCT's Note

$B$ が下に凸な時、それを以下のように並べた行列はMongeである。

B = { 9, 4, 1, 0, 1, 4 }
↓
i\j  0  1  2  3  4  5  6 ...
 0   9 ∞ ∞ ∞ ∞ ∞ ∞      これはMonge
 1   4  9 ∞ ∞ ∞ ∞ ∞
 2   1  4  9 ∞ ∞ ∞ ∞
 3   0  1  4  9 ∞ ∞ ∞
 4   1  0  1  4  9 ∞ ∞
 5   4  1  0  1  4  9 ∞
 6  ∞  4  1  0  1  4  9
 7  ∞ ∞  4  1  0  1  4
 :

で、Mongeな行列に対し、列全体に同じ値を加算してもMonge性は保たれるという性質がある。

A = { 3, 1, 4, 1, 5, 9, 2 }

i\j  0  1  2  3  4  5  6
    +3 +1 +4 +1 +5 +9 +2
 0  12 ∞ ∞ ∞ ∞ ∞ ∞      これもMonge
 1   7 10 ∞ ∞ ∞ ∞ ∞
 2   4  5 13 ∞ ∞ ∞ ∞
 3   3  2  8 10 ∞ ∞ ∞
 4   4  1  5  5 14 ∞ ∞
 5   7  2  4  2  9 18 ∞
 6  ∞  5  5  1  5 13 11
 7  ∞ ∞  8  2  4 10  6

これは、$X_{i,j}=A_j+B_{i-j}$ の行列になっている。これで行毎に最小値を取ったものは、まさしく $C$ である。

$B$ のみが下に凸な場合、畳み込みの添字に注意。$X_{i,j}=A_{i-j}+B_j$ としてしまうと上手くいかない。

• LARSCH Algorithm | Library
• 簡易版 LARSCH Algorithm - noshi91のメモ

SMAWKは、あらかじめ全ての $A$ が（関数で部分的に求めるにしろ）わかっていないといけない。

しかし動的計画法で順番に求める必要がある場合に、それだと上手くいかないことがある。

LARSCHは、「$\arg\min_k A_{j,k}$ の値が求まって、はじめて $A_{i,j+1}$ の値を使用できる」場合でも求めることができる。