Mergeselect

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(n2). 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 x0 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 s0 udává kolik prvků je menších než mnou zvolený „pivot“ x0. Tady nastane rekurze třetího kroku. Pokud je s0 větší než m, zvolím x1 menší než x0, pokud je s0 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) * log2(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(n2) (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(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/diskrimi­nace, 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í/para­lelizovatelné 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í.

Flattr this!

This entry was posted in Algo, řazení. Bookmark the permalink.

5 Responses to Mergeselect

  1. v6ak says:

    „Chci najít m-tý nejmenší element, jinými slovy potřebuji najít prvek, kterému předchází m menších bratří.“

    Není to spíše prvek, kterému předchází m-1 menších bratří?

    Nemám teď úplně pocit, že jsem to pochopil, ale:
    • Co když k=1? Bude i zde platit tato slozitost?
    • Pokud je asymptotická slozitost založená na pravděpodobnosti, nedostaneš worst-case, ale spíše průměr. V tomto směru na tom ale je stejně jako QuickSelect, ne?

    • m nebo m-1, záleží jestli indexy začínají od 0 nebo 1. Ale máš pravdu, m-1 dává smysl.

      k=1 je stejně patologický případ jako k=n a pro ně to fungovat nebude.

      Přesně tak, nejde o nejhorší případ, ale o ten očekávaný, stejně jako quickselect. Doplnil jsem to do textu.

  2. v6ak says:

    Pokud je řeč o počtu prvků, pak je jedno, jestli začínají od nuly, nebo od jednicky.

    No, k ani nemůže být zvoleno tak, aby se obecně rovnalo n, pokud není počet prvků předem znám. A pokud je, pak vše děláme v O(1), až na výjimky, kde máme O(∞). Pak ale asymptotická složitost vlastně říká jen to, zda je garantováno nebo časově omezeno dokončení algoritmu. Pokud bychom k vybrali az podle délky vstupu, uz by to nebyla konstanta, takže by slo o jiný algoritmus.

    Varianta k=1 by ale nemusela být patologická, pokud sekvence nebude začínat největším prvkem. (Podobný problém může nastat i u většího k, jen s nižší pravděpodobností.)

    • Když k = n, pak to běží v čase n log(n), protože udělám jeden průchod mergesortu, seřadím vstupní kolekci a hotovo vyberu m-tý prvek.

      Když k = 1, stane se z toho svým způsobem quickselect. První blok velikosti 1 je pivot, použitý pro filtrování ostatních singulárních bloků.

      • v6ak says:

        k = n (nebo obecneji k≥n): Jak se to vezme. Pokud se pred spustenim podivam na delku kolekce, a takto zvolim k, pak ano, ale je to uz jiny algoritmus. Pokud jde o to, ze pro nekdere male vstupy to takto vyjde, tak to nema na asymptotickou analyzu zadny vliv :)

        k=1: Asi ano, ale z hlediska asymptoticke slozitosti to patologicke neni, ne? Prumer bude porad O(n) a worst-case je IMHO dosazitelny i u vetsich k, pokud se ty nejvetsi prvku dostanou do prvniho bloku.

Leave a Reply

Your email address will not be published. Required fields are marked *