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

Quick­se­lect všichni známe a mi­lu­jeme – jde o al­go­rit­mus pro na­le­zení k-tého nejmen­šího prvku ne­se­řa­ze­ného pole v prů­měr­ném čase O(n). Jediný jeho ne­do­sta­tek je, že v pa­to­lo­gic­kých pří­pa­dech může běžet v čase O(n2). Quicksort to má ale stejně a je to sou­část dohody, kterou jsme od­sou­hla­sili, když jsme pus­tili do quick-věcí. Pro ty, kteří chtějí ga­rance, exis­tuje o něco ra­fi­no­va­nější al­go­rit­mus medián me­di­ánů, který má ga­ran­to­vaný O(n) čas a ne­pří­liš pří­z­nivé kon­stantní fak­tory.

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

Když jsem před ně­ja­kou dobou při­pra­vo­val před­nášku o řazení, trávil jsem víc času, než by bylo zá­hodno, du­má­ním nad pa­ra­le­lami, po­dob­nostmi a du­a­lismy. Nedá moc práce a člověk si začne všímat ana­lo­gií mezi al­go­ritmy, často se lišící jen pře­vrá­ce­nými šip­kami. Bádání mě o pár týdnů poz­ději do­vedlo k na­črt­nutí níže uve­de­ného di­a­gramu2 . Každá jeho osa před­sta­vuje nějaký prin­cip, způsob fun­go­vání nebo atri­but. V místě, kde se osy pro­tí­nají, exis­tuje al­go­rit­mus, ne­mu­sím ho znát, ale vždycky tam je.

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

A tady se do­stá­vám k ti­tulku článku. Pokud sy­me­t­rie platí a exis­tuje zr­ca­dlový vztah mezi quicksor­tem a merge sortem, musí také exis­to­vat merge obdoba quick­se­lectu – tedy hy­po­te­tický mer­ge­se­lect, který stejně jako jeho quick-bratr běží v li­ne­ár­ním čase.

Ně­ja­kou dobu jsem nad touhle otáz­kou – za­jí­ma­vou, i když z prak­tic­kého hle­diska ne­vý­znam­nou – pře­mýš­lel, hnán touhou po uspo­řá­da­ném vesmíru, který dává smysl a ale nemohl jsem přijít na řešení. Ale pak jed­noho osa­mě­lého od­po­ledne jsem zažil heu­réka moment. O(n) mer­ge­se­lect fun­guje, pokud se do­pus­tím všech hříchů asympto­tické ana­lýzy, budu před­stí­rat, že kon­stanty ne­hrají žádnou roli1 , a v jednom místě taktně přejdu ná­strahy prav­dě­po­dob­nosti.

Re­kur­zivní al­go­rit­mus mer­ge­se­lectu pra­cuje ve čtyřech kro­cích zná­zor­ně­ných na ná­sle­du­jí­cím sché­matu.

V prvním kroku roz­dě­lím vstup na bloky kon­stantní délky k > 1 (např. 1024), kte­rých je tedy n/k.

Ve druhém začnu s merge sortem, ne­do­táhnu ho však do konce, ale sor­tuji jen uvnitř bloků. Vý­sled­kem jsou se­řa­zené bloky délky k, kdy každý z nich je zpra­co­ván v čase k * log(k). Slo­ži­tost tohoto kroku je tedy: (n / k) * (k * log(k)) = n * log(k), ale pro­tože k je kon­stanta, od­po­vídá to O(n).

Třetí krok je místo, kde se děje všechna magie. Chci najít m-tý nejmenší ele­ment, jinými slovy po­tře­buji najít prvek, kte­rému před­chází m-1 men­ších bratří. Začnu tím, že si zvolím pro­střední ele­ment x0 prv­ního se­řa­ze­ného bloku a bi­nár­ním hle­dání ve všech ostat­ních blo­cích najdu stejný nebo nej­bližší větší prvek. V prů­běhu pro­hle­dá­vání sčítám pozice těchto prvků. Součet s0 udává kolik prvků je men­ších než mnou zvo­lený „pivot“ x0. Tady na­stane re­kurze tře­tí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. Pro­vá­dím tedy další bi­nární hle­dání.3

Slo­ži­tost tře­tího kroku je (n / k) * log(k) * log(k) = (n / k) * log2(k). Stále platí, že k je kon­stanta, a výraz tedy od­po­vídá O(n).

Pokud platí, že sx = m, našel jsem m-tý nejmenší prvek. To však nebude příliš prav­dě­po­dobné, pro­tože se hle­daný ele­ment bude s nej­větší prav­dě­po­dob­ností na­chá­zet v ostat­ních blo­cích.

Ve čtvr­tém kroku se na­chá­zím v si­tu­aci, kdy jsem ve všech blo­cích bi­nár­ním hle­dá­ním určil horní a dolní mez, kde by se m-tý ele­ment měl na­chá­zet a tyto dílčí sub­sek­vence jsou se­řa­zené. Co udělám? Začnu je mer­go­vat do bloků ma­xi­mální délky k a celý al­go­rit­mus opa­kuji, ne­hle­dám však m-tý nejmenší prvek, ale ten na pozici (m-s), kde s je součet dol­ních mezí. Pokud je vý­sled­kem jen jeden blok, můžu v něm přímo najít kýžený prvek.

Aby vý­sledná re­kurze měla li­ne­ární čas, je třeba splnit je­di­nou pod­mínku – množ­ství práce v další ite­raci se musí zmen­šit o kon­stantní ná­so­bek 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 po­skládá do O(n).

Jaká je velké f v pří­padě mer­ge­se­lectu? Tady přijde na řadu zmi­ňo­vaný úkrok stra­nou kolem prav­dě­po­dob­nosti. Sub-sek­vence vy­me­zená dolní a horní mezí v prvním bloku má délku 0 nebo 1 (buď ob­sa­huje kýžený prvek nebo ne). Když budu uva­žo­vat uni­formní roz­dě­lení dat, ve­li­kosti sub-sek­vencí ve všech ostat­ních blo­cích by měly být po­dobně velké jako to v prvním bloku.4 Mělo by tedy dojít k re­dukci cel­kové ve­li­kosti na něco jako n/k. Pokud tohle platí, pak je oče­ká­vaná slo­ži­tost celé re­kurze mer­ge­se­lectu na úrovni O(n), nej­horší případ má slo­ži­tost O(n2) (stejně jako quick­se­lect) a to je celkem fajn.

Do­datky:

Pseu­do­kó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:

Po­známky:

  1. Pře­kva­pivě hrají. Quicksort byl pu­b­li­ko­ván roku 1961 a nej­vět­ších po­kroků bylo do­sa­ženo zlep­šo­vá­ním právě kon­stant­ních fak­torů. Pů­vodní Ho­a­rova metoda zvolit jako pivot první (nebo ná­hod­nou) po­ložku má prů­měrný kon­stantní faktor 1.386, metoda nej­lepší ze tří 1.188 a kvazi-nej­lepší-z-devíti 1.094. Další pokrok byl three way quicksort, který je odolný proti pa­to­lo­gic­kým vstu­pům.
  2. S nej­větší prav­dě­po­dob­nostní na něm chybí jedna až tři osy. Tři hlavní osy na di­a­gramu jsou merge/dis­kri­mi­nace, po­rov­nání/radix, dělení na dva/široké. Další osy by měly být ještě efemérní proces/dy­na­mická datová struk­tura, sek­venční/pa­ra­le­li­zo­va­telné a možná i líné/striktní. Ale když začnu za­chá­zet to ta­ko­vých de­tailů, hra­nice se roz­ma­zá­vají a víc než tech­nický obor to začne při­po­mí­nat te­o­lo­gické ar­gu­menty.
  3. Pro­tože k je při­ja­telně velká kon­stanta je možné bi­nární hle­dání pro­vá­dět velice efek­tivně, jak je po­psáno v Binary Search eli­mi­na­tes Branch Mispre­dicti­ons.
  4. A to, podle ně­ko­lika si­mu­lací, které jsem udělal, vypadá, že platí.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články