Mergeselect

22. 8. 2016

Quickselect všichni známe a milujeme - jde o algoritmus pro nalezení k-tého nejmenšího prvku neseřazeného pole v průměrném čase O(n). Jediný jeho nedostatek je, že v patologických případech může běžet v čase O(n²). Quicksort to má ale stejně a je to součást dohody, kterou jsme odsouhlasili, když jsme pustili do quick-věcí. Pro ty, kteří chtějí garance, existuje o něco rafinovanější algoritmus medián mediánů, který má garantovaný O(n) čas a nepříliš příznivé konstantní faktory.

To všechno víme už 55 let.

Když jsem před nějakou dobou připravoval přednášku o řazení, trávil jsem víc času, než by bylo záhodno, dumáním nad paralelami, podobnostmi a dualismy. Nedá moc práce a člověk si začne všímat analogií mezi algoritmy, často se lišící jen převrácenými šipkami. Bádání mě o pár týdnů později dovedlo k načrtnutí níže uvedeného diagramu² . Každá jeho osa představuje nějaký princip, způsob fungování nebo atribut. V místě, kde se osy protínají, existuje algoritmus, nemusím ho znát, ale vždycky tam je.

Quicksort a merge sort jsou dva nečekaní bratři, ležící na opačných koncích jedné osy.

A tady se dostávám k titulku článku. Pokud symetrie platí a existuje zrcadlový vztah mezi quicksortem a merge sortem, musí také existovat merge obdoba quickselectu - tedy hypotetický mergeselect, který stejně jako jeho quick-bratr běží v lineárním čase.

Nějakou dobu jsem nad touhle otázkou - zajímavou, i když z praktického hlediska nevýznamnou - přemýšlel, hnán touhou po uspořádaném vesmíru, který dává smysl a ale nemohl jsem přijít na řešení. Ale pak jednoho osamělého odpoledne jsem zažil heuréka moment. O(n) mergeselect funguje, pokud se dopustím všech hříchů asymptotické analýzy, budu předstírat, že konstanty nehrají žádnou roli¹ , a v jednom místě taktně přejdu nástrahy pravděpodobnosti.

Rekurzivní algoritmus mergeselectu pracuje ve čtyřech krocích znázorněných na následujícím schématu.

V prvním kroku rozdělím vstup na bloky konstantní délky k > 1 (např. 1024), kterých je tedy n/k.

Ve druhém začnu s merge sortem, nedotáhnu ho však do konce, ale sortuji jen uvnitř bloků. Výsledkem jsou seřazené bloky délky k, kdy každý z nich je zpracován v čase k * log(k). Složitost tohoto kroku je tedy: (n / k) * (k * log(k)) = n * log(k), ale protože k je konstanta, odpovídá to O(n).

Třetí krok je místo, kde se děje všechna magie. Chci najít m-tý nejmenší element, jinými slovy potřebuji najít prvek, kterému předchází m-1 menších bratří. Začnu tím, že si zvolím prostřední element x₀ prvního seřazeného bloku a binárním hledání ve všech ostatních blocích najdu stejný nebo nejbližší větší prvek. V průběhu prohledávání sčítám pozice těchto prvků. Součet s₀ udává kolik prvků je menších než mnou zvolený "pivot" x₀. Tady nastane rekurze třetího kroku. Pokud je s₀ větší než m, zvolím x₁ menší než x₀, pokud je s₀ větší než m, zvolím naopak. Provádím tedy další binární hledání.³

Složitost třetího kroku je (n / k) * log(k) * log(k) = (n / k) * log²(k). Stále platí, že k je konstanta, a výraz tedy odpovídá O(n).

Pokud platí, že s_x = m, našel jsem m-tý nejmenší prvek. To však nebude příliš pravděpodobné, protože se hledaný element bude s největší pravděpodobností nacházet v ostatních blocích.

Ve čtvrtém kroku se nacházím v situaci, kdy jsem ve všech blocích binárním hledáním určil horní a dolní mez, kde by se m-tý element měl nacházet a tyto dílčí subsekvence jsou seřazené. Co udělám? Začnu je mergovat do bloků maximální délky k a celý algoritmus opakuji, nehledám však m-tý nejmenší prvek, ale ten na pozici (m-s), kde s je součet dolních mezí. Pokud je výsledkem jen jeden blok, můžu v něm přímo najít kýžený prvek.

Aby výsledná rekurze měla lineární čas, je třeba splnit jedinou podmínku - množství práce v další iteraci se musí zmenšit o konstantní násobek f, tedy n_i+1 = f * n_i. f nemusí být nijak veliké, stačí když je o něco málo větší než 1 a výsledný čas se poskládá do O(n).

Jaká je velké f v případě mergeselectu? Tady přijde na řadu zmiňovaný úkrok stranou kolem pravděpodobnosti. Sub-sekvence vymezená dolní a horní mezí v prvním bloku má délku 0 nebo 1 (buď obsahuje kýžený prvek nebo ne). Když budu uvažovat uniformní rozdělení dat, velikosti sub-sekvencí ve všech ostatních blocích by měly být podobně velké jako to v prvním bloku.⁴ Mělo by tedy dojít k redukci celkové velikosti na něco jako n/k. Pokud tohle platí, pak je očekávaná složitost celé rekurze mergeselectu na úrovni O(n), nejhorší případ má složitost O(n²) (stejně jako quickselect) a to je celkem fajn.

Dodatky:

Tenhle algoritmus jsem nikde neviděl. To ale moc neznamená, protože jsem jednak diletant, který nerozumí asymptotické analýze, a druhak jsem nehledal nijak důkladně. Kdybych si myslel, že jsem objevil něco zcela nového, tak by to nepříjemně zavánělo naivitou a arogancí. Už párkrát jsem byl přesvědčený, že jsem na stopě něčeho inovativního, jen abych o týden později našel v paperu z roku 1975 zmínku o té mojí novince. Navíc nešlo o stěžejní myšlenku publikace, ale jen drobnou noticku ve stylu: "Jo a tohle je tady věc, o které se ví."
Mergeselect (stejně jako standardní merge sort) není in-place. Potřebuje jednu horní a jednu dolní mez pro každý blok, kterých je n/k. To odpovídá O(n) extra prostoru, což je víc než O(log(n)), který je považován za in-place. Parametr k může však být relativně velký a proto můžou být praktické nároky snesitelné.
Algoritmus nenaprogramoval, protože nechci ztratit pár dnů debugováním krajních případů.
Merge sort je možné nahradit Timsortem nebo dokonce zcela jiným O(n log(n)) řadícím algoritmem, takže vztah mergeselectu k merge sortu je spíš jen duchovní.
Po vzoru quickselectu se dá zkonstruovat radixselect, který místo podle pivotu dělí dle radixu.

Pseudokód:

def mergeselect[T](xs: Seq[T], m: Int): T = {
  require(m < xs.length)

  // #1
  val k = 1024
  val blocks = xs.grouped(k)

  // #2
  for (block <- blocks) {
    mergesort(block)
  }

  if (blocks.length == 1) {
    return blocks(m)
  }

  // #3
  val lowerBounds = Array.fill(blocks.length)(0)
  val upperBounds = Array.tabulate(blocks.length) { b => b.length }
  val positions = Array.fill(blocks.length)(0)

  // binary search for s == m
  var s = 0
  while (lowerBounds(0) <= upperBounds(0)) {
    val p = (lowerBounds(0) + upperBounds(0)) >>> 1 // mid pisition
    val x = blocks(0)(p) // mid value

    s = p

    for (b <- 1 until blocks.length) {
      val block = blocks(b)
      val pos = binarySearch(x, in = block, from = lowerBounds(b), to = upperBounds(b))
      positions(b) = pos
      s += pos
    }

    if (s < m) {
      lowerBounds(0) = p + 1
    } else if (s > m) {
      upperBounds(0) = p - 1
    } else {
      return x
    }

    if (s < m) {
      for (b <- 1 until blocks.length) {
        lowerBounds(b) = positions(b) + 1
      }
    } else {
      for (b <- 1 until blocks.length) {
        upperBounds(b) = positions(b) - 1
      }
    }
  }

  // #4
  val newXs: Seq[T] = 0 until blocks.length map { b => blocks(b).slice(from = lowerBounds(b), to = upperBounds(b)) }.flatten
  mergeselect(newXs, m-s)
}

Dále k tématu:

Poznámky:

Překvapivě hrají. Quicksort byl publikován roku 1961 a největších pokroků bylo dosaženo zlepšováním právě konstantních faktorů. Původní Hoarova metoda zvolit jako pivot první (nebo náhodnou) položku má průměrný konstantní faktor 1.386, metoda nejlepší ze tří 1.188 a kvazi-nejlepší-z-devíti 1.094. Další pokrok byl three way quicksort, který je odolný proti patologickým vstupům.
S největší pravděpodobnostní na něm chybí jedna až tři osy. Tři hlavní osy na diagramu jsou merge/diskriminace, porovnání/radix, dělení na dva/široké. Další osy by měly být ještě efemérní proces/dynamická datová struktura, sekvenční/paralelizovatelné a možná i líné/striktní. Ale když začnu zacházet to takových detailů, hranice se rozmazávají a víc než technický obor to začne připomínat teologické argumenty.
Protože k je přijatelně velká konstanta je možné binární hledání provádět velice efektivně, jak je popsáno v Binary Search eliminates Branch Mispredictions.
A to, podle několika simulací, které jsem udělal, vypadá, že platí.