Вывод оптимального алгоритма с помощью формализма Бёрда-Меертенса

Некоторые оптимальные алгоритмы, оказывается, можно вывести из неоптимальных, пользуясь эквивалентными преобразованиями алгоритма. Бёрд и Меертенс разработали формализм, который устанавливает свойства функций высшего порядка map, fold, scan, позволяющие преобразовывать алгоритмы в эквивалентные. (См. также на Вики). Ниже представлен вольный перевод статьи Бёрда.

Рассмотрим задачу поиска максимальной суммы сегмента массива. Эту задачу можно переформулировать в виде математически точного ответа:

Для всех сегментов, которые можно получить из массива, необходимо посчитать сумму чисел, а затем среди всех таких сумм найти максимальную.

Такая словесная формулировка неудобна для последующих манипуляций. Поэтому введём дополнительные определения.

Функция segs возвращает все возможные сегменты (подмассивы) массива:

def segs(input: Array[Int]): List[Array[Int]] = ???

segs :: [Int] -> [[Int]]

Функция sum вычисляет сумму элементов:

def sum(input: Array[Int]): Int = input.sum

sum :: [Int] -> Int

Чтобы получить все возможные суммы, нам достаточно применить к каждому подмассиву функцию sum. Такая операция делается функцией map:

def map[A, B](f: A => B)(list: List[A]): List[B] = list.map(f)

map :: (a -> b) -> [a] -> [b]

Теперь, чтобы найти максимальную сумму, достаточно применить функцию max:

def max(lst: List[Int]): Int = lst.max

max :: [Int] -> Int

Обозначим функцию, находящую решение задачи, mss ("max segment sum"). Эта функция будет иметь примерно такой вид:

def mss(arr: Array[Int]): Int = max(map(sum)(segs(arr)))

mss :: [Int] -> Int
mss arr = max(map(sum)(segs(arr)))

В Haskell'е можно ту же самую функцию записать элегантнее, если использовать оператор композиции .. Такой оператор возвращает функцию одного аргумента, которая вначале применит правую функцию к своему аргументу, а затем левую функцию к результату, возвращённому правой функцией. Иными словами, f . g эквивалентно x => f(g(x)).

mss = max . map(sum) . segs

Также можно опустить скобочки, т.к. аргументы передаются и так:

mss = max . map sum . segs

В Scala вместо . используется compose:

val mss: Array[Int] => Int = 
    max `compose` map(sum) `compose` segs

Полученное выражение является одной из возможных формулировок исходной задачи и в то же время даёт алгоритм её решения. Этот алгоритм имеет сложность O(n^3) и требует O(n^2) памяти.

Вслед за Бёрдом попробуем произвести эквивалентные преобразования этого алгоритма.

Определение segs

Получить все сегменты можно в два приёма. Вначале получить все возможные начала массива (inits). А затем для каждого начала найти все возможные окончания (tails).

def inits(arr: Array[Int]): List[Array[Int]] = ???
def tails(arr: Array[Int]): List[Array[Int]] = ???

inits :: [Int] -> [[Int]]
tails :: [Int] -> [[Int]]

Чтобы применить tails к каждому началу, воспользуемся функцией map:

segs1 :: [Int] -> [[[Int]]]
segs1 arr = map(tails)(inits(arr))
segs1 = map tails . inits

К сожалению, результатом будет вложенный список списков сегментов. Чтобы получить простой список, склеим все вложенные списки с помощью функции concat:

concat :: [[a]] -> [a]

def concat[A](lst2: List[List[A]]): List[A] = lst2.flatten

Таким образом, функция segs может быть реализована так:

segs = concat . map tails . inits

val segs = concat `compose` map(tails) `compose` inits

Подставив это определение в исходный алгоритм, получим:

mss = max . map sum . concat . map tails . inits

val mss: Array[Int] => Int = 
    max `compose` map(sum) `compose` concat `compose` map(tails) `compose` inits

Распределение map внутри вложенных списков (map promotion)

Если мы склеиваем вложенные списки, а затем каждый элемент преобразуем с помощью какой-то функции, то мы с тем же успехом можем вначале преобразовать каждый элемент, а уже затем склеивать вложенные списки. Результат будет тем же.

Чтобы преобразовывать каждый элемент внутри вложенных списков, нам потребуется дважды использовать map. Один раз для внешнего списка, а потом, для каждого вложенного.

def mapconcat[A, B](f: A => B)(lst2: List[List[A]]): List[B] = 
  lst2.map(innerList => innerList.map(f)).flatten
def mapconcat2[A, B](f: A => B)(lst2: List[List[A]]): List[B] = 
  lst2.flatten.map(f)

mapconcat f lst2 = concat(map(map(f))(lst2))
mapconcat f = concat . map(map(f))
mapconcat f = concat . map (map f)

mapconcat2 f lst2 = map(f) (concat(lst2))
mapconcat2 f = map(f) . concat

Эти две функции дают одинаковый результат. То есть они эквивалентны. Будем обозначать такого рода эквивалентность символом ≡:

map(f) . concat ≡ concat . map(map f) 

Воспользуемся этим свойством и преобразуем mss:

mss = max . map sum . concat . map tails . inits
mss = max . concat . map (map sum) . map tails . inits

val mss: Array[Int] => Int = 
    max `compose` map(sum) `compose` concat `compose` map(tails) `compose` inits
val mss: Array[Int] => Int = 
    max `compose` concat `compose` map(map(sum)) `compose` map(tails) `compose` inits

Распределение max по вложенным спискам с использованием распределительного свойства concat/fold (fold promotion)

Функция max, используемая выше, не является элементарной. Она принимает на вход список элементов, которые потом обрабатывает последовательно, на каждом шагу применяя бинарную операцию максимизации. Такую операцию можно обозначить стрелочкой ⋏ (в дальнейшем мы будем использовать термин "операция max", чтобы было понятнее). Функция max, принимающая на вход список элементов, может быть реализована с помощью бинарной операции max и функции свёртки foldl:

max = foldl(⋏)(-∞)

def max(lst: List[Int]): Int =
  lst.foldLeft(Int.MinValue)(math.max)

Операция max обладает полезным свойством ассоциативности. Можно вначале посчитать max для вложенных списков, а потом уже посчитать max для результатов. Иными словами:

max . concat ≡ max . map(max)

(Здесь мы сразу используем символ эквивалентности ≡.)

def maxConcat(lst2: List[List[Int]]): Int = 
  max(concat(lst2))
def maxConcat2(lst2: List[List[Int]]): Int = 
  max(lst2.map(innerList => max(innerList)))

Такое же свойство работает для любых ассоциативных операций и позволяет избавиться от склейки списков.

Воспользуемся этим свойством и преобразуем mss:

mss = max . concat  . map (map sum) . map tails . inits
mss = max . map max . map (map sum) . map tails . inits

val mss: Array[Int] => Int = 
    max `compose` concat `compose` map(map(sum)) `compose` map(tails) `compose` inits
val mss: Array[Int] => Int = 
    max `compose` map(max) `compose` map(map(sum)) `compose` map(tails) `compose` inits

Дистрибутивный закон для map и compose (слияние циклов)

Если мы последовательно применяем функции к элементам списка, то это можно сделать либо двумя проходами map которые склеиваются compose, либо одним проходом, применяя склеенную функцию:

map f . map g ≡ map (f . g)

В последнем выражении для mss у нас имеется три подряд map'а. Можно склеить все три:

map max . map (map sum) . map tails ≡ map (max . (map sum) . tails)

Таким образом, получаем новую версию mss:

mss = max . map max . map (map sum) . map tails . inits
mss = max . map (max . (map sum) . tails)       . inits

val mss: Array[Int] => Int = 
    max `compose` map(max) `compose` map(map(sum)) `compose` map(tails) `compose` inits
val mss: Array[Int] => Int = 
    max `compose` map(    max `compose` map(sum) `compose` tails      ) `compose` inits

Схема Горнера

До сих пор правила преобразований были не очень сложными и затрагивали отдельные операции. Схема Горнера позволяет упростить свёртку с двумя операциями сразу.

В схеме Горнера задействованы две операции ⊕, ⊗. Для операции ⊗ необходим нейтральный элемент, обозначаемый 1, кроме того, должен выполняться дистрибутивный закон относительно ⊕. Запишем вначале эту схему для этих двух операций:

Мы видим, что в исходном выражении имеется n*(n-1)/2 операций ⊗, а в последнем — только n. Т.е. такая схема выгодна с вычислительной точки зрения.

Применим эту схему к нашему алгоритму. В качестве операции ⊗ у нас будет обычное сложение (+), нейтральным элементом 1 для которого является 0. А в качестве операции ⊕ мы будем использовать ⋏ (операцию max). Для сложения выполняется дистрибутивный закон относительно операции max:

(То есть, можно найти максимум, а затем прибавить константу, либо вначале прибавить контстанту, а потом найти максимум.) Для этих операций схема Горнера принимает вид:

Нам удобнее перенумеровать переменные в обратном порядке:

Запишем левую часть с использованием функций, рассмотренных выше. В каждой скобке находятся элементы массива, начиная с некоторого индекса и до конца. Все элементы для всех скобок предоставляются функцией tails.

В скобках производится суммирование элементов. Чтобы произвести суммирование в каждом наборе элементов, воспользуемся функцией map с параметром sum:

map sum . tails

И все полученные скобки затем сворачиваются с использованием операции max. Т.е. необходимо применить функцию max:

max . map sum . tails

Рассмотрим теперь правую часть схемы Горнера. Здесь перебираются все элементы слева направо и постепенно вмешиваются в аккумулятор с помощью пары операций — сложения и max. Таким же образом работает foldl (foldLeft). Нам надо лишь выделить операцию, которая будет выполняться каждый раз (UPD: спасибо CorwinH, обратившему внимание на то, что аккумулятор идёт первым аргументом):

(accumulator, element) => element + accumulator ⋏ 0

Обозначим такую операцию ⊙ (или summax0). Тогда

val summax0 =
  (accumulator: Int, element: Int) => element + math.max(accumulator, 0)

⊙ accumulator element = element + (max 0 accumulator)
⊙ = (+) . (max 0)

Начало вычислений следует начинать с нейтрального элемента для операции max. В качестве такого нейтрального элемента возьмём -∞ (или Int.MinValue). Правая часть схемы Горнера принимает вид:

foldLeft(Int.MinValue)(summax0)

foldl(⊙)(-∞)

И в целом схема Горнера:

max . map sum . tails ≡ foldl(⊙)(-∞)

Воспользуемся схемой, чтобы упростить выражение для mss:

mss = max . map (max . (map sum) . tails) . inits
mss = max . map (foldl(⊙)(-∞)           ) . inits

val mss: Array[Int] => Int = 
    max `compose` map(    max `compose` map(sum) `compose` tails      ) `compose` inits
val mss: Array[Int] => Int = 
    max `compose` map(    foldLeft(Int.MinValue)(summax0)             ) `compose` inits

Сканирование с накоплением

Функции foldl и foldr осуществляют поэлементную свёртку с аккумулятором и возвращают самое последнее значение. Иногда требуется вернуть все промежуточные результаты, получаемые в процессе свёртки, в виде списка. Такие функции также имеются в стандартной библиотеке и носят названия scanl и scanr (scanLeft, scanRight):

scanl :: (b -> a -> b) -> b -> [a] -> [b]

Обычная функция foldl легко выражается через scanl. Достаточно взять последний элемент:

foldl f zero ≡ last . scanl f zero

Также можно и scanl выразить через foldl. Достаточно просто вызвать foldl для всех начал списка:

scanl f zero ≡ map (foldl f zero) . inits

Эти два эквивалентных алгоритма отличаются производительностью. В правой части сложность составляет O(n^2), а в левой — O(n). Поэтому это соотношение имеет смысл применять в обратную сторону, для улучшения эффективности. В частности, используя в качестве функции f операцию summax0 (⊙), а в качестве начального значения Int.MinValue (-∞), получим такое соотношение:

map (foldl ⊙ -∞) . inits
≡
scanl ⊙ -∞

map(foldLeft(Int.MinValue)(summax0)) `compose` inits
≡
scanLeft(Int.MinValue)(summax0)

Подставляя в выражение для mss, получаем:

mss = max . map (foldl(⊙)(-∞)) . inits
mss = max . scanl ⊙ -∞

val mss: Array[Int] => Int = 
    max `compose` map(foldLeft(Int.MinValue)(summax0)) `compose` inits
val mss: Array[Int] => Int = 
    max `compose` scanLeft(Int.MinValue)(summax0)

Свёртка сканирования

В нашем последнем выражении мы дважды проходим по списку длиной n и по дороге вынуждены сохранять промежуточные результаты в памяти O(n). Можем совместить эти два прохода и сохранять промежуточные результаты для каждого из них в двух отдельных ячейках. Для этого достаточно в качестве аккумулятора использовать пару.

Запишем функцию max через foldl:

max ≡ foldl ⋏ -∞

Аккумулятор для этой части алгоритма будем хранить в первом элементе пары.

scanl возвращает значения своего аккумулятора на каждом шаге в виде списка. Т.к. мы формируем один проход, то нам больше не потребуется список, а только лишь значение аккумулятора на каждом шаге. Этот аккумулятор будем хранить во втором элементе пары.
Сконструируем функцию maxsummax0, которая будет обрабатывать пару аккумуляторов и следующий элемент (UPD по замечанию CorwinH):

def maxsummax0(pairOfAccumulators: (Int, Int), element: Int): (Int, Int) =
  val (w1, v1) = pairOfAccumulators
  val v2 = summax0(v1, element)
  val w2 = math.max(w1, v2)
  (w2, v2)

maxsummax0 (w1, v1) element = (w2, v2)
  where
    v2 = summax0 v1 element
    w2 = max w1 v2

Для foldl также требуется начальное значение аккумуляторов. Т.к. начальное значение scanl равно 0 и мы обязательно применяем к нему операцию max (⋏), то начальным значением для foldl будет пара нулей (0, 0).

В конце алгоритма нам нужен только результат максимизации. Т.е. первый элемент пары fst (_._1).

Таким образом, получаем следующее соотношение

max . scanl ⊙ -∞ ≡ fst . foldl maxsummax0 (0,0)

def fst[A, B](p: (A, B)): A = p._1

fst . foldLeft((0,0))(maxsummax0)

Используем это соотношение в mss:

mss = max . scanl ⊙ -∞
mss = fst . foldl maxsummax0 (0,0)

val mss: Array[Int] => Int = 
    max `compose` scanLeft(Int.MinValue)(summax0)
val mss: Array[Int] => Int = 
    fst . foldLeft((0,0))(maxsummax0)

Заключение

Начав с буквальной формулировки постановки задачи, которая по совместительству оказалась алгоритмом, имеющим сложность O(n^3), Бёрд в результате ряда эквивалентных преобразований получил линейный алгоритм (O(n)):

mss = max . map sum . segs
≡
mss = fst . foldl maxsummax0 (0,0)

Получившийся алгоритм является оптимальным для этой задачи алгоритмом Кадане. Вывод алгоритма заодно является математическим доказательством того, что этот алгоритм является решением поставленной задачи.