funkcionálně.cz

Přední český blog o funkcionálním programování, kde se o funkcionálním programování nepíše
««« »»»

Mergeselect

22. 8. 2016 — k47

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^2). 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 diagramu2 . 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 roli1 , 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_0 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_0 udává kolik prvků je menších než mnou zvolený "pivot" x_0. Tady nastane rekurze třetího kroku. Pokud je s_0 větší než m, zvolím x_1 menší než x_0, pokud je s_0 větší než m, zvolím naopak. Provádím tedy další binární hledání.3

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

Pokud platí, že sx = 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 ni+1 = f ni. 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.4 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^2) (stejně jako quickselect) a to je celkem fajn.

Dodatky:

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(for = 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:

  1. 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.
  2. 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.
  3. 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.
  4. A to, podle několika simulací, které jsem udělal, vypadá, že platí.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články