funkcionálně.cz

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

Jak řadit v lineárním čase, křísit mrtvé a dosáhnout osvícení

11. 4. 2016

Když se řekne řazení, většina lidí, kteří napsali víc než jednu smyčku a dva ify se zeptá: "A je to nutné?" Řazení nemá pověst nejrychlejší operace pod sluncem a je preferováno, když se tomu dá vyhnout např. hashováním. Někdy je to ale nutné a v tom případ je třeba zatnout zuby a připravit se na nějaký ten rozpálený křemík.

Jedinci, kterým na zdi visí portrét Donalda Knutha si okamžitě vybaví, že řazení má v nejlepším případě časovou složitost O(n log(n)) a vzpomenou si na všechny ty quicksorty, merge sorty, heap sorty a různé monstrózní hybridy.

Pravda je ale o něco komplikovanější: Řadící algoritmy založené na porovnání potřebují provést řádově n * log(n) porovnání. Existují však i jiné způsoby řazení, které jsou založené na frekvenci, rozložení a histogramech, a ty nepotřebují porovnávat nic s ničím. Stačí jim lineární čas, a co víc: Jsou zatraceně rychlé.

Jedním z nejjednodušších je counting sort, který funguje tak, že vstupní pole jednou projdu, zaznamenávám kolikrát jsem na který element narazil a na konci položky vyplivnu ve požadovaném pořadí.

Může to vypadat například takhle (psáno v jazyce fast-and-loose Scala):

def countingSort(arrayToBeSorted: Array[Int]): Array[Int] = {
  val buckets = new ArrayBuffer[Int]()

  for (el <- arrayToBeSorted) {
    buckets(el) += 1 // increment counter
  }

  val result = new ArrayBuffer[Int]()

  for {
    i <- 0 until buckets.size // element
    _ <- 0 until buckets(i)   // count
  } result += i // append element

  result.toArray
}

Bystří si jistě už všimli, že slibovaná lineární složitost má pár nedostatků. Jak je vidět, je třeba mít po ruce pole, které je velké jako největší element ve vstupní kolekci. Co se stane, když bych se pokusil seřadit čísla 1 a 18446744073709551616? Musel bych přikoupit hodně paměti.

Counting sort má lineární časovou složitost, ale potřebuje prostor přímo úměrný rozpětí možných hodnot. To je ideální, když řadím například pole bajtů, ale v ostatních případech to má k ideálu velice daleko.

Řešením je radix sort. Ten funguje velice podobně, ale k řazení nepoužívá celou hodnotu najednou (tj. všechny dostupné bity), ale vždy jen její část. Tak například, když bych chtěl seřadit pole čtyřbajtových intů, radix sort může v první iteraci data rozdělit podle nejdůležitějšího bajtu. Vznikne tak 256 skupin (partition), kdy každá z nich obsahuje jen elementy shodující se v nejvyšším bajtu. V další iteraci provedu stejný krok pro každou skupinu zvlášť s druhým nejdůležitějším bajtem. Vznikne tak 256x256 podskupin, z nichž každá sdílí horních 16 bitů. Když udělám 4 rekurzivní iterace a spojím vzniklé podoblasti, dostanu seřazenou sekvenci. Stačí mi k tomu jen 4 iterace vstupním polem a pokud se nepletu, 4 lze považovat za konstantu.

Časová složitost radix sortu je O(kn), kde n je délka vstupu a k je velikost jedné položky. Když řadím pole primitivních typů nebo malých structů, k je malé a všechno běží velice rychle. Když se však rozhodnu radixovat dlouhé stringy nebo velké a komplikované (nedej bože nekonečné) objekty, k začne dominovat a situace je o poznaní méně růžová. Zvlášť s přihlédnutím k faktu, že log2(232), což je maximální rozsah 32 bitového intu, je jen 32, k tedy nemá moc místa pro manévrování. Stačí když radixovaný objekt trochu nakyne a najednou n log(n) představuje správnou cestu vpřed.

Majíce tohle všechno na paměti jsem napsal vlastní implementaci radix sortu ve Scale, abych na vlastní oči viděl jak rychlé je rychlé řazení. Můj kód je založený na moudrosti tohoto článku a jde o LSD radix sort. Neřadí tedy od nejvýznamnějšího bajtu, ale od toho nejméně významného. To se může zdát neintuitivní, ale funguje to, protože radix sort je stabilní algoritmus - více detailů ve zmiňovaném článku. Pro čtyřbajtové inty potřebuje provést 5 iterací nad daty - v jedné spočítá histogram frekvencí jednotlivých bajtů a ve čtyřech zbývajících přesouvá položky podle histogramů.

To však není jediný způsob. Existuje spousta dalších odrůd radix sortu: je možné řadit od nejdůležitějších bitů, od těch nejméně důležitých, může být in-place, může potřebovat extra O(n) místa, může být hybridní (použije se jiný algoritmus, když je jasné, že se nevyplatí radixovat), může být paralelní, může brát v potaz cache, může používat software managed buffers, atd. atd.

Výsledek? Jak řekl klasik: ho-ly-shit. Přesvědčte se sami:

Graf ukazuje rychlost řazení polí obsahujících náhodné inty mým radix sortem a funkcí java.util.Arrays.sort, což je hybrid mezi dual-pivot quicksortem, mergesortem a insertion sortem. Jak je vidět, radix sort běží v lineárním čase a navíc je velice rychlý. A to jde o neaktuální graf, který měří rychlost verze bez několika optimalizací3 , a navíc do průměru započítává i počáteční pomalé běhy, kdy JVM nebyla zahřátá. Malé zpomalení s rostoucí velikostí dat není způsobeno asymptotami algoritmu, ale hardwarovými efekty. Všechny operace jsou rychlé, čtení je vždycky sekvenční, zápis je skoro sekvenční do 256 míst v paměti, které se pohodlně vejdou do cache5 , neobsahuje žádné nepředvídatelné skoky, které by vedly k brach-miss - situace je zcela odlišná od quicksortu, který neobsahuje nic jiného než nepředvídatelné skoky.1

1GB dat (250M zcela náhodných intů) je na mém stroji možné seřadit za 4 vteřiny, a 1GB floatů za 2.9 vteřiny2 . To odpovídá 14.9 ns, resp 10.8 ns na každý element pole. Když se celé pole vejde do cache, řazení je o něco rychlejší. Každému intu i floatu stačí 9.7 ns, což odpovídá rychlosti řazení přes 100M/s v jednom vlákně. Radix sort navíc umí poznat, když není třeba řadit podle určitého bajtu a přeskočí celý jeden řadící krok, což vede ke značnému zrychlení. Například když řadím pole obsahující jen čísla mezi 0 a 64k, dvě iterace z pěti jsou přeskočeny a je třeba jen 8 ns místo 14.9 ns na element.

Řazení floatů radix sortem je poměrně delikátní záležitost, která se spoléhá na to, že čísla mají bitovou reprezentaci podle specifikace IEEE754. Nejjednodušší možností je data seřadit jako integery a pak opravit pořadí záporných hodnot. O něco rafinovanější je floaty na vstupu transformovat, tak aby by je bylo možné řadit jako 32-bitová čísla bez znaménka a na konci je změnit zpátky na float. Nejlepší řešení je však tomu předejít a upravit offsety záporných radixů.

Ok, takže když je radix sort tak rychlý, k čemu se dá použít? K řazení, pochopitelně. Ale najde i uplatnění tam, kde se řazení zdá, jako jít s kanónem na vrabce.

Jedno takové místo je seskupování.

Když bych chtěl seskupit moře hodnot podle jejich klíče z páru klíč-hodnota (kde klíč i hodnota jsou čtyřbajtové inty zabalené do osmibajtového longu), mohl bych to s do očí bijící neefektivitou napsat v jednom řádku Scaly.

keyvals.groupBy(getKey).map { case (k, vs) => vs.map(getKey) }

To je pěkné, ale nešlo by to rychleji, nejlépe s minimem alokací? Šlo. S pomocí mapy, která obsahuje měnitelné buffery.

val buffers = mutable.Map[Int, mutable.ArrayBuilder.ofInt]()

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  buffers.getOrElseUpdate(k, new mutable.ArrayBuilder.ofInt) += v
}

buffers.toSeq.sortBy(_._1).values

Když zanedbám režii mapy (která může být zanedbatelná, ale také nemusí), za pozornost stojí poslední řádek: sortBy. Najednou jsem tam, kde jsem začal. Když je velikost každé skupiny malá, je jejich počet úměrný velikosti vstupních dat a konečné řazení bude velice drané.

Když vím, že klíče náleží do určitého rozpětí, můžou použít pole bufferů.

val buffers = new Array[mutable.ArrayBuilder.ofInt](length)

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  if (buffers(k) == null) {
    buffers(k) = new mutable.ArrayBuilder.ofInt
  }

  buffers(k) += v
}

buffers.filter(_ != null).map(_.result)

Jak je vidět, na konci není třeba nic řadit, protože výsledek je přirozeně seřazen podle klíče ve stylu counting sortu. Pokud mě trápí opakované alokace a kopírování dat v rostoucích bufferech a můžu si dovolit dva průchody daty, je možné v jenom průchodu spočítat velikosti jednotlivých skupin, alokovat pro ně pole velká tak jak je potřeba, a v druhém průchodu je naplnit. To může být v některých případech rychlejší a v jiných, kdy nechci přijít s OutOfMemoryError po funuse, nezbytné.

Tohle je lepší, ale nešlo by to... však víte... rychleji?

Tady do hry vstupuje radix sort. Co kdybych pole seřadil podle prvních 4 bajtů (tj. té části, která obsahuje klíč) a pak výsledek v jedné iteraci rozřezal na kousíčky. Nějak takhle:

val scratch = new Array[Long](keyvals.length)
val (sorted, _) = RadixSort.sort(keyvals, scratch, 0, keyvals.length, 4, 8, false)

var i = 0
(0 until numberOfKeys iterator) map { key =>
  var start = i

  // find end of the current group
  while (i < keyvals.length && getKey(sorted(i)) == key) { i += 1 }
  val end = i

  // empty group
  if (start == end) {
    null

  // copy values  into a separate array
  } else {
    val res = new Array[Int](end - start)
    var j = 0
    while (start < end) {
      res(j) = getVal(sorted(start))
      start += 1
      j += 1
    }

    (key, res)
  }
} filter (_ != null )

Na rozdíl od předchozích případů je kód o něco složitější, na druhou stranu nejde o pseudokód, ale skutečně fungující fragment. A co víc: Je skutečně rychlejší než všechny předchozí ukázky.

Proč?

Jedno vysvětlení je nasnadě: Pokud je skupin hodně a data jsou více méně náhodná, v předchozím případě každý nový klíč může vést ke několika drahým cache-miss. Jeden načte více méně náhodný pointer na buffer a druhý načte objekt bufferu a třetí konec pole v bufferu na které se bude vkládat. A už jsem zmiňoval, že jít do DRAM je drahé? Myslím, že ano. Radix sort tyhle problémy nemá, párkrát projede data, seřadí je a výsledek se nakonec ještě jednou proletí a rozseká na kousíčky. A všechno udělá velice rychle.

Dále k tématu:

  1. Špatně odhadnuté skoky představují tak velký problém, že někdy je quicksort rychlejší, pokud je záměrně zvolen špatný pivot. Když pivot data rozdělí fifty-fifty, je skok nepředvídatelný. Když je ale rozdělí třeba 25-75, procesor bude odhadovat, že skok půjde doprava a většinou se trefí. viz: An Experimental Study of Sorting and Branch Prediction. V tomto kontextu za pozornost také stojí: Branch Prediction and the Performance of Interpreters - Don't Trust Folklore a Branchless code sequences.
  2. Předpokládám, že v případě floatů JVM JIT zkompiluje kód s využitím SIMD instrukcí a proto je výsledek rychlejší. Nebo také může jít o korelace mezi bitovými vzory v reprezentaci floatů. Podrobnosti jsem nezkoumal.
  3. Jedna z optimalizací, která může za velice specifických situací vést k dvojnásobnému zrychlení je zpětná iterace jedné smyčky. V běžných situacích, kdy jsou vstupní data o něco větší než cache, přidá k dobru ~10%. Když je vstupní pole veliké, nemá žádný vliv. Další optimalizaci, kterou jsem zkoušel, byly software managed buffers zmiňované zde a zde. Ty fungují tak, že data nezapisuji přímo do paměti, ale do malého bufferu, který se vejde do cache a do paměti ho zkopíruji jen, když se naplní. Tohle přidalo pár procent výkonu pro řazení intů, ale vedlo to k drastickému zpomalení v případě floatů. Výsledek nebyl nijak přesvědčivý a proto jsem tuhle optimalizaci nakonec nepoužil.
  4. To s přihládnutím k tomu, že radix sort musí udělat 5 průchodů celým polem znamená 2.2 ns nebo 7 taktů CPU na element v každé iteraci. Z paměti je přitom třeba číst/zapisovat 4.8 GB/s (vysvětlení v tomto paperu, který hlásí, že jeden průchod jejich radix sortu, který dává 240M/s na čtyřech jádrech, trvá 13.2 taktu na element - pozor jejich čísla nesedí).
  5. V literatuře se uvádí, že radix sort nedělá sekvenční zápis v tom smyslu, že nezapisuje do jednoho místa, ale do mnoha, v tomto případě 256. To vede k neoptimálnímu využití propustnosti pamětí a může vést k neustálým výpadkům TLB. Moje experimenty na moderním hardwaru však nejsou přesvědčivé aspoň, co se týká TLB. Perf hlásí skoro nulu dTLB-load-misses a dTLB-store-misses, ale na druhou stranu cache miss se stále nějaké vyskytují (jako load-miss tak i store-miss).
  6. Se složitostí radix sortu je to o něco komplikovanější. Je pravda, že běží v čase O(kn), ale povaha onoho k je trochu matoucí. Když bych měl n různých čísel, nejmenší počet bitů potřebných k reprezentaci každého z nich je log<sub>2</sub>(n). Ono O(kn) je tedy vlastně jen O(n log(n)) v přestrojení.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články