秦九韶アルゴリズムを用いて最大部分配列和問題を解く

23/12/2023

最大部分配列和問題
秦九韶のアルゴリズム
Kadaneのアルゴリズム
元の推論
可愛さの裏側

WARNING

この記事はChatGPTによって中国語から翻訳されたもので、いくつかの誤りが含まれているかもしれません。不正確な部分があればご了承ください。

問有兩尖田一段，其尖長不等，兩大斜三十九步，兩小斜二十五步，中廣三十步，欲知其積幾何？
術曰：以少廣求之，翻法入之。置半廣自乘為半冪，與小斜冪相減相乘為小率，以半冪與大斜冪相減相乘為大率。以二率相減餘自乘為實，併二率倍之為從上廉，以一為益隅，開翻法三乘方得積。

部の人々は、関数型プログラミングを軽視しています。彼らは、関数型プログラミングが無意味な抽象化を増やすだけで、その代償はパフォーマンスの低下だと感じるかもしれません。このような見解を持つ人々は、おそらくRichard Birdの著書 PEARLS OF FUNCTIONAL ALGORITHM DESIGN を読んだことがないでしょう。この本の内容は確かに難解で、私も一度も完全に読み通したことはありません。

幸運なことに、Richard Birdは1989年にわずか5ページの論文 Algebraic identities for program calculation を書いています。この論文では、最大部分配列和¹（Maximum subarray problem）という有名な問題を用いて、関数型プログラミング思考がアルゴリズム問題に対する独特の理解を見事に示しています：最大部分配列和問題は秦九韶アルゴリズムで解くことができます。

❖最大部分配列和問題

最大部分配列和問題の目的は、数列の一次元方向で連続した部分配列を見つけ、その部分配列の和が最大になるようにすることです。例えば、数列 [-2, 1, -3, 4, -1, 2, 1, -5, 4] に対して、連続した部分配列の中で最大の和を持つのは [4, -1, 2, 1] で、その和は 6 です²。

明らかに、私たちは平凡なアルゴリズムを持っています、すなわち、すべての部分配列を列挙し、最大値を求めます：

1def mss(arr):
2    r = 0
3    for end in range(len(arr) + 1):
4        for start in range(0, end):
5            s = 0
6            for k in range(start, end):
7                s += arr[k]
8            r = max(s, r)

このプログラムはHaskellで書くとシンプルになります：

1import Data.List (tails, inits)
2
3mss = maximum . map sum . segs
4segs = concat . map tails . inits

ここで tails は「接尾辞」、inits は「接頭辞」を意味します：

1tails [1, 2, 3] = [[1, 2, 3], [2, 3], [3], []]
2inits [1, 2, 3] = [[], [1], [1, 2], [1, 2, 3]]

明らかに、「すべての接頭辞のすべての接尾辞」、または「すべての接尾辞のすべての接頭辞」はリストのすべての部分リストになります。

さて、アルゴリズム自体に戻りましょう。この単純なアルゴリズムの時間複雑度は $O(n^3)$ です。すべての部分配列を列挙するのに $O(n^2)$ 、各部分配列の和を求めるのに $O(n)$ が必要なので、結果として $O(n^3)$ になります。

実際には、最も優れたアルゴリズムは $O(n)$ の時間複雑度でこの問題を解決できます。Birdの論文は、単純なアルゴリズムから出発して、最も優れたアルゴリズムを_導き出す_方法について議論しています。まず、アルゴリズムを少し変形します：

1def mss(arr):
2    m = 0
3    for end in range(len(arr) + 1):
4        m = max(max_tails(arr, end), m)
5
6def max_tails(arr, end):
7    m = 0
8    for start in range(0, end):
9        s = 0
10        for k in range(start, end):
11            s += arr[k]
12        m = max(m, s)
13    return m

これは単純な変形のように見えますが、私は3つのループの内側の2つを取り出しただけです。max_tails 関数は、配列 arr の [0, end) の範囲の最大の接尾辞の和を返します。最大値関数を $\uparrow$ で表すと、max_tails の関数バージョン、または数学的バージョンは次のように書けます：

\begin{aligned} \textsf{arr[0:end]} &\equiv [ x_1, x_2, \cdots, x_n ] \\ m &= (\sum_{i=1}^n x_i) \uparrow (\sum_{i=2}^n x_i) \uparrow \cdots \uparrow (x_n) \uparrow 0 \end{aligned}

これは一見すると非常に難解に見えますが、以下の例で説明します。まず、max_tails([1, 2, 3], 3)、つまり [1, 2, 3] のすべての接尾辞の和の最大値を考えてみましょう。

max_tails の計算は次のように書けます：

1(1 + 2 + 3) ↑ (2 + 3) ↑ (3) ↑ 0

つまり、max_tails は各接尾辞の和を計算し、それらの最大値を求めます。最大値を二項関数 $\uparrow$ として書くと、それを使って各接尾辞の和をつなげて最終的な式を構成することが自然になります。

朴素アルゴリズムでは、接尾辞の最大値を求めるために $O(n^2)$ の複雑さが必要です。以下では、秦九韶のアルゴリズムを用いてこの複雑さを $O(n)$ に減らします。

❖秦九韶のアルゴリズム

秦九韶は彼の著作「数術九章」で、高次方程式の近似解を求める方法を提案しました。その一部は現代の研究者によって「秦九韶のアルゴリズム」と呼ばれています。

$n$ 次多項式を考えてみましょう。

f(x) = \sum_{i=0}^{n} a_i x^{i}

特定の $t$ に対して、 $f(t)$ の値を求める方法はいくつかあります。最も単純な方法は、各項の値を個別に求めてから合計することです。このアルゴリズムでは $O(n^2)$ 回の乗算が必要です。しかし、秦九韶のアルゴリズムは次のように注意しています：

\begin{aligned} & a_n x^n + a_{n-1} x^{n-1} + \cdots + a_0 \\ = & ((a_n x^{n - 1}) + (a_{n-1} x^{n - 2}) + \cdots + a_{1}) x + a_0 \\ = & \cdots \\ = & (((a_n x + a_{n - 1}) x + \cdots) x + a_{1})x + a_0 \end{aligned}

これは反復形式で書くことができます：

\begin{aligned} s_0 &= a_n \\ s_{i + 1} &= s_i t + a_{n - i - 1} \end{aligned}

このように、 $s_n = f(t)$ となります。この計算では、乗算は $n$ 回、加算も $n$ 回行われます。

一見すると、秦九韶のアルゴリズムは単なる 数値計算 のアルゴリズムに過ぎません。それは max_tails 関数と何の関係があるのでしょうか？ここでは、「一般化」された秦九韶のアルゴリズム形式を考える必要があります。

まず、このアルゴリズムの入力は多項式である必要はありません。以下の式は秦九韶のアルゴリズムで書き換えることができます：

\begin{aligned} & x_1x_2x_3 + x_2x_3 + x_3 + 1 \\ =& (((x_1 + 1) x_2) + 1) x_3 + 1 \end{aligned}

次に、秦九韶のアルゴリズムが考慮する演算は $(+, \cdot)$ である必要はありません。このアルゴリズムが公因数を抽出するとき、実際には乗法分配法則を用いて式を書き換えています。任意の二つの関数 $(\oplus, \odot)$ を考えると、それらが以下の条件³を満たす限り、秦九韶のアルゴリズムを適用することができます：

$\odot$ は $\oplus$ に対して（右）分配します
$\begin{aligned} a \cdot c + b \cdot c &= (a + b) \cdot c \\ (a \odot c) \oplus (b \odot c) &= (a \oplus b) \odot c \\ \end{aligned}$
$\odot$ の（左）単位元が存在します、つまり、すべての $x$ に対して $a \odot x = x$ となる $a$ が存在します。

以下の式を例にとります：

\begin{aligned} & ((x_1 \odot x_2) \odot x_3) \oplus (x_2 \odot x_3) \oplus (x_3) \oplus (1) \\ = & ((((x_1 \oplus 1) \odot x_2) \oplus 1) \odot x_3) \oplus 1 \end{aligned}

ここで $1$ は $\odot$ の左単位元です。

❖Kadaneのアルゴリズム

先ほど述べたように、max_tails は次のように書くことができます：

\begin{aligned} \textsf{arr[0:end]} &\equiv [ x_1, x_2, \cdots, x_n ] \\ m &= (\sum_{i=1}^n x_i) \uparrow (\sum_{i=2}^n x_i) \uparrow \cdots \uparrow (x_n) \uparrow 0 \end{aligned}

ここでの $\uparrow$ は $\oplus$ と見なすことができ、 $+$ は $\odot$ と見なすことができます。加法には単位元 $0$ が存在し、分配性も明らかです：

(a \uparrow b) + c = (a + c) \uparrow (b + c)

続けて [1, 2, 3] を例にとります，

1(1 + 2 + 3) ↑ (2 + 3) ↑ (3) ↑ 0
2= ((1 + 2) ↑ 2 ↑ 0) + 3 ↑ 0
3= (((1 ↑ 0) + 2) ↑ 0) + 3 ↑ 0
4= ((((0 + 1) ↑ 0) + 2 ↑ 0) + 3) ↑ 0

この性質に基づいて、秦九韶のアルゴリズムを用いて max_tails を書き換えます：

1def step(s, a):
2    # (s + a) ↑ 0
3    return max(s + a, 0)
4
5def max_tails(arr, end):
6    s = 0
7    for i in range(end):
8        s = step(s, arr[i])

明らかに、書き換えた後の max_tails の時間複雑度は O(n) です。さらに、次のことに注意します：

\begin{aligned} & \textsf{max\_tails}(\textsf{arr}, \textsf{end} + 1) \\ =\ & \textsf{step}(\textsf{step}(\cdots,\textsf{arr}[\textsf{end} - 1]), \textsf{arr}[\textsf{end}]) \\ =\ & \textsf{step}(\textsf{max\_tails(arr, end)}, \textsf{arr}[\textsf{end}]) \end{aligned}

mssが [1, end + 1) の max_tails を求める必要があるため、改写後の時間複雑度は $O(n^2)$ となります。しかし、上記の発見から、max_tails(arr, i) の値は max_tails(arr, i - 1) の値から求めることができます。これをすぐに $O(n)$ の mss 関数を得るために使用することができます：

1def step1(state, a):
2    (s, m) = state
3    s_next = step(s, a)
4    return (s_next, max(m, s_next))
5
6def mss(arr):
7    state = (0, 0)
8    for a in arr:
9        state = step1(state, a)
10    return state[1]

間違いなく、これは大名鼎鼎の Kadane アルゴリズムです。

❖元の推論

Bird 教授は彼の研究を一つの段落にまとめています：

私は、明確だが非効率的な関数型プログラム、つまり手元の問題の仕様として機能するプログラムを取り、等式推論を用いてより効率的なものを計算するという特定のタスクに興味がありました。⁴

つまり、Bird 教授は、シンプルで明確だが効率的でない関数型プログラムから出発し、一般的なルールを用いてプログラムを最適化し、より効率的なプログラムを得ることを好んでいます。Bird 教授は最も直感的なプログラムから出発し、一歩一歩推論を進めることで、KMP アルゴリズムのような有名なアルゴリズムを得ることができます。Bird 教授の論文では、最大部分配列和問題に対する元の推論は次のようになっています⁵：

1mss = { by definition }                                     
2      max . map sum . segs                                  O(n³)
3    = { by definition of segs }
4      max . map sum . concat . map tails . inits            O(n³)
5    = { map promotion }
6      max . concat . map (map sum) . map tails . inits      O(n³)
7    = { definition of max and fold promotion }               
8      max . map max . map (map sum) . map tails  . inits    O(n³)
9    = { map distributivity }
10      max . map (max . (map sum) . tails) . inits           O(n³)
11    = { Horner's Rule }
12      max . map (foldl (⊗) 0) . inits                       O(n²)
13    = { scan theorem }
14      max . scanl (⊗) 0                                     O(n)
15    = { fold-scan fusion }
16      fst . foldl (⊙) (0, 0)                                O(n)
17    where a ⊗ b = (a + b) ↑ 0
18          (u, v) ⊙ t = (w ↑ u, w), w = v ⊗ t

たった8つの等式で最終的なアルゴリズムを得ることができ、その技術は驚異的と言えます。さらに価値あることは、秦九韶のアルゴリズムや最終的な最適化が、推論過程で直接的な定理として示されていることです。これらの定理は、関数型プログラミングでよく使われる関数（map、foldl、concatなど）の等式に関するものです。例えば、秦九韶のアルゴリズムは次のように直接表現することができます：

1fold (⊕) b . map (fold (⊗) a) . tails = fold (⊙) a
2u ⊙ v = (u ⊗ v) ⊕ a

では、なぜ関数型プログラミングはこのようなプログラムの推論を容易にするのでしょうか？それは、関数型プログラミングが私たちにプログラムの 代数的 性質により注意を向けさせるからです。Bird 教授はまた、The algebra of programmingという本も書いています。

❖可愛さの裏側

読者が上記のコードを LeetCode に提出しようとすると、いくつかのテストで失敗するかもしれません。これはもちろん、アルゴリズムが間違っているわけではなく、最大部分配列和の問題の表現が異なるからです。

(LeetCode): 配列が与えられたとき、その最大部分配列和を出力します。
(Programming Pearls): 配列が与えられたとき、その最大部分配列和を出力します。もし配列の最大部分配列和が負であれば、0を返します。

私たちの問題は Programming Pearls のバージョンで、最大部分配列和が負になることはありません。

LeetCode バージョンの問題を処理するには、以下のような簡単な修正が可能です：

1mssNeg l
2  | r == 0    = max l
3  | otherwise = r
4    where r = mss l

この方法は少し重たい感じがします。Bird 教授の導出過程を修正して、負数のケースを処理できるアルゴリズムを導出することは可能でしょうか？この問題に対して、私は理想的な解決策を提供していません。なぜなら、Bird 教授の最初の式が間違っているからです。

1--sum [] = 0
2mss = maximum . map sum . segs

segs の結果には必ず空のリストが存在し、その空のリストを sum で評価すると、0が得られます。したがって、この式は負の最大和のケースを処理することはできません。また、ここでの秦九韶のアルゴリズムでは、0 は + の単位元であるため、負数を処理できるアルゴリズムに修正するには、完全に一からやり直す必要があるかもしれません。

これは、関数型プログラムの導出の可愛らしさの裏側かもしれません。関数型プログラムの導出はプログラムの代数的性質に依存していますが、この代数的性質は時として脆弱で、問題を少し修正するだけで、元の代数的性質が消失することがあります。命令型プログラムにとっては、それ自体がこれらの代数的性質に依存していないため、些細な変更がすぐに機能するかもしれませんが、関数型、特に関数型プログラムの導出から導出されたプログラムにとっては、このような変更は些細なものではありません。

Bird 教授にメールで意見を聞いてみてもいいですか？

それもいいかもしれません。しかし、Bird 教授は2022年4月4日に永遠に私たちを去りました。天国に手紙を送ることができる郵便局があるのでしょうか……