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 — k47

Když se řekne řazení, vět­šina lidí, kteří na­psali víc než jednu smyčku a dva ify se zeptá: „A je to nutné?“ Řazení nemá pověst nej­rych­lejší ope­race pod slun­cem a je pre­fe­ro­váno, když se tomu dá vy­hnout např. ha­sho­vá­ním. Někdy je to ale nutné a v tom případ je třeba za­tnout zuby a při­pra­vit se na nějaký ten roz­pá­lený křemík.

Je­dinci, kterým na zdi visí por­trét Do­nalda Knutha si oka­mžitě vybaví, že řazení má v nej­lep­ším pří­padě ča­so­vou slo­ži­tost O(n log(n)) a vzpo­me­nou si na všechny ty quicksorty, merge sorty, heap sortyrůzné mon­strózní hyb­ridy.

Pravda je ale o něco kom­pli­ko­va­nější: Řadící al­go­ritmy za­lo­žené na po­rov­nání po­tře­bují pro­vést řádově n * log(n) po­rov­nání. Exis­tují však i jiné způ­soby řazení, které jsou za­lo­žené na frek­venci, roz­lo­žení a his­to­gra­mech, a ty ne­po­tře­bují po­rov­ná­vat nic s ničím. Stačí jim li­ne­ární čas, a co víc: Jsou za­tra­ceně rychlé.

Jedním z nej­jed­no­duš­ších je coun­ting sort, který fun­guje tak, že vstupní pole jednou projdu, za­zna­me­ná­vám ko­li­krát jsem na který ele­ment na­ra­zil a na konci po­ložky vy­plivnu ve po­ža­do­va­ném pořadí.

Může to vy­pa­dat na­pří­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 sli­bo­vaná li­ne­ární slo­ži­tost má pár ne­do­statků. Jak je vidět, je třeba mít po ruce pole, které je velké jako nej­větší ele­ment ve vstupní ko­lekci. Co se stane, když bych se po­ku­sil se­řa­dit čísla 1 a 18446744073709551616? Musel bych při­kou­pit hodně paměti.

Coun­ting sort má li­ne­ární ča­so­vou slo­ži­tost, ale po­tře­buje pro­stor přímo úměrný rozpětí mož­ných hodnot. To je ide­ální, když řadím na­pří­klad pole bajtů, ale v ostat­ních pří­pa­dech to má k ideálu velice daleko.

Ře­še­ním je radix sort. Ten fun­guje velice po­dobně, ale k řazení ne­po­u­žívá celou hod­notu na­jed­nou (tj. všechny do­stupné bity), ale vždy jen její část. Tak na­pří­klad, když bych chtěl se­řa­dit pole čtyř­baj­to­vých intů, radix sort může v první ite­raci data roz­dě­lit podle nej­dů­le­ži­těj­šího bajtu. Vznikne tak 256 skupin (par­ti­tion), kdy každá z nich ob­sa­huje jen ele­menty sho­du­jící se v nej­vyš­ším bajtu. V další ite­raci pro­vedu stejný krok pro každou sku­pinu zvlášť s druhým nej­dů­le­ži­těj­ším bajtem. Vznikne tak 256x256 pod­sku­pin, z nichž každá sdílí hor­ních 16 bitů. Když udělám 4 re­kur­zivní ite­race a spojím vzniklé podob­lasti, do­stanu se­řa­ze­nou sek­venci. Stačí mi k tomu jen 4 ite­race vstup­ním polem a pokud se ne­pletu, 4 lze po­va­žo­vat za kon­stantu.

Časová slo­ži­tost radix sortu je O(kn), kde n je délka vstupu a k je ve­li­kost jedné po­ložky. Když řadím pole pri­mi­tiv­ních typů nebo malých structů, k je malé a všechno běží velice rychle. Když se však roz­hodnu ra­di­xo­vat dlouhé stringy nebo velké a kom­pli­ko­vané (nedej bože ne­ko­nečné) ob­jekty, k začne do­mi­no­vat a si­tu­ace je o po­znaní méně růžová. Zvlášť s při­hléd­nu­tím k faktu, že log2(232), což je ma­xi­mální rozsah 32 bi­to­vého intu, je jen 32, k tedy nemá moc místa pro ma­né­v­ro­vání. Stačí když ra­di­xo­vaný objekt trochu nakyne a na­jed­nou n log(n) před­sta­vuje správ­nou cestu vpřed.


Majíce tohle všechno na paměti jsem napsal vlastní im­ple­men­taci radix sortu ve Scale, abych na vlastní oči viděl jak rychlé je rychlé řazení. Můj kód je za­lo­žený na moud­rosti tohoto článku a jde o LSD radix sort. Neřadí tedy od nej­vý­znam­něj­šího bajtu, ale od toho nejméně vý­znam­ného. To se může zdát ne­in­tu­i­tivní, ale fun­guje to, pro­tože radix sort je sta­bilní al­go­rit­mus – více de­tailů ve zmi­ňo­va­ném článku. Pro čtyř­baj­tové inty po­tře­buje pro­vést 5 ite­rací nad daty – v jedné spo­čítá his­to­gram frek­vencí jed­not­li­vých bajtů a ve čtyřech zbý­va­jí­cích pře­souvá po­ložky podle his­to­gramů.

To však není jediný způsob. Exis­tuje spousta dal­ších odrůd radix sortu: je možné řadit od nej­dů­le­ži­těj­ších bitů, od těch nejméně dů­le­ži­tých, může být in-place, může po­tře­bo­vat extra O(n) místa, může být hyb­ridní (po­u­žije se jiný al­go­rit­mus, když je jasné, že se ne­vy­platí ra­di­xo­vat), může být pa­ra­lelní, může brát v potaz cache, může po­u­ží­vat soft­ware ma­naged bu­f­fers, atd. atd.

Vý­sle­dek? Jak řekl klasik: ho-ly-shit. Pře­svědčte se sami:

Graf uka­zuje rych­lost řazení polí ob­sa­hu­jí­cích ná­hodné inty mým radix sortem a funkcí java.util.Arrays.sort, což je hybrid mezi dual-pivot quicksor­tem, mer­gesor­tem a in­ser­tion sortem. Jak je vidět, radix sort běží v li­ne­ár­ním čase a navíc je velice rychlý. A to jde o ne­ak­tu­ální graf, který měří rych­lost verze bez ně­ko­lika op­ti­ma­li­zací3 , a navíc do prů­měru za­po­čí­tává i po­čá­teční pomalé běhy, kdy JVM nebyla za­hřátá. Malé zpo­ma­lení s ros­toucí ve­li­kostí dat není způ­so­beno asympto­tami al­go­ritmu, ale hard­wa­ro­vými efekty. Všechny ope­race jsou rychlé, čtení je vždycky sek­venční, zápis je skoro sek­venční do 256 míst v paměti, které se po­ho­dlně vejdou do cache5 , ne­ob­sa­huje žádné ne­před­ví­da­telné skoky, které by vedly k brach-miss – si­tu­ace je zcela od­lišná od quicksortu, který ne­ob­sa­huje nic jiného než ne­před­ví­da­telné skoky.1

1GB dat (250M zcela ná­hod­ných intů) je na mém stroji možné se­řa­dit za 4 vte­řiny, a 1GB floatů za 2.9 vte­řiny2 . To od­po­vídá 14.9 ns, resp 10.8 ns na každý ele­ment pole. Když se celé pole vejde do cache, řazení je o něco rych­lejší. Kaž­dému intu i floatu stačí 9.7 ns, což od­po­vídá rych­losti řazení přes 100M/s v jednom vlákně. Radix sort navíc umí poznat, když není třeba řadit podle ur­či­tého bajtu a pře­skočí celý jeden řadící krok, což vede ke znač­nému zrych­lení. Na­pří­klad když řadím pole ob­sa­hu­jící jen čísla mezi 0 a 64k, dvě ite­race z pěti jsou pře­sko­čeny a je třeba jen 8 ns místo 14.9 ns na ele­ment.

Řazení floatů radix sortem je po­měrně de­li­kátní zá­le­ži­tost, která se spo­léhá na to, že čísla mají bi­to­vou re­pre­zen­taci podle spe­ci­fi­kace IEEE754. Nej­jed­no­dušší mož­ností je data se­řa­dit jako in­te­gery a pak opra­vit pořadí zá­por­ných hodnot. O něco ra­fi­no­va­nější je floaty na vstupu trans­for­mo­vat, tak aby by je bylo možné řadit jako 32-bitová čísla bez zna­ménka a na konci je změnit zpátky na float. Nej­lepší řešení je však tomu pře­de­jít a upra­vit off­sety zá­por­ných radixů.


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

Jedno takové místo je se­sku­po­vání.

Když bych chtěl se­sku­pit moře hodnot podle jejich klíče z páru klíč-hod­nota (kde klíč i hod­nota jsou čtyř­baj­tové inty za­ba­lené do os­mi­baj­to­vého longu), mohl bych to s do očí bijící ne­e­fek­ti­vi­tou napsat v jednom řádku Scaly.

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

To je pěkné, ale nešlo by to rych­leji, nej­lépe s mi­ni­mem alo­kací? Šlo. S pomocí mapy, která ob­sa­huje mě­ni­telné bu­f­fery.

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ž za­ne­dbám režii mapy (která může být za­ne­dba­telná, ale také nemusí), za po­zor­nost stojí po­slední řádek: sortBy. Na­jed­nou jsem tam, kde jsem začal. Když je ve­li­kost každé sku­piny malá, je jejich počet úměrný ve­li­kosti vstup­ních dat a ko­nečné řazení bude velice drané.

Když vím, že klíče náleží do ur­či­tého rozpětí, můžou použít pole bu­f­ferů.

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, pro­tože vý­sle­dek je při­ro­zeně se­řa­zen podle klíče ve stylu coun­ting sortu. Pokud mě trápí opa­ko­vané alo­kace a ko­pí­ro­vání dat v ros­tou­cích bu­f­fe­rech a můžu si do­vo­lit dva prů­chody daty, je možné v jenom prů­chodu spo­čí­tat ve­li­kosti jed­not­li­vých skupin, alo­ko­vat pro ně pole velká tak jak je po­třeba, a v druhém prů­chodu je na­pl­nit. To může být v ně­kte­rých pří­pa­dech rych­lejší a v jiných, kdy nechci přijít s Ou­tO­f­Me­mo­ryError po funuse, ne­zbytné.

Tohle je lepší, ale nešlo by to… však víte… rych­leji?

Tady do hry vstu­puje radix sort. Co kdy­bych pole se­řa­dil podle prv­ních 4 bajtů (tj. té části, která ob­sa­huje klíč) a pak vý­sle­dek v jedné ite­raci roz­ře­zal na kou­síč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řed­cho­zích pří­padů je kód o něco slo­ži­tější, na druhou stranu nejde o pseu­do­kód, ale sku­tečně fun­gu­jící frag­ment. A co víc: Je sku­tečně rych­lejší než všechny před­chozí ukázky.

Proč?

Jedno vy­svět­lení je na­snadě: Pokud je skupin hodně a data jsou více méně ná­hodná, v před­cho­zím pří­padě každý nový klíč může vést ke ně­ko­lika drahým cache-miss. Jeden načte více méně ná­hodný poin­ter na buffer a druhý načte objekt bu­f­feru a třetí konec pole v bu­f­feru na které se bude vklá­dat. A už jsem zmi­ňo­val, že jít do DRAM je drahé? Myslím, že ano. Radix sort tyhle pro­blémy nemá, pár­krát pro­jede data, seřadí je a vý­sle­dek se na­ko­nec ještě jednou pro­letí a roz­seká na kou­síčky. A všechno udělá velice rychle.


Dále k tématu:


  1. Špatně od­had­nuté skoky před­sta­vují tak velký pro­blém, že někdy je quicksort rych­lejší, pokud je zá­měrně zvolen špatný pivot. Když pivot data roz­dělí fifty-fifty, je skok ne­před­ví­da­telný. Když je ale roz­dělí třeba 25-75, pro­ce­sor bude od­ha­do­vat, že skok půjde do­prava a vět­ši­nou se trefí. viz: An Ex­pe­ri­men­tal Study of Sor­ting and Branch Pre­diction. V tomto kon­textu za po­zor­nost také stojí: Branch Pre­diction and the Per­for­mance of In­ter­pre­ters – Don't Trust FolkloreBranchless code sequen­ces.
  2. Před­po­klá­dám, že v pří­padě floatů JVM JIT zkom­pi­luje kód s vy­u­ži­tím SIMD in­strukcí a proto je vý­sle­dek rych­lejší. Nebo také může jít o ko­re­lace mezi bi­to­vými vzory v re­pre­zen­taci floatů. Po­drob­nosti jsem ne­zkou­mal.
  3. Jedna z op­ti­ma­li­zací, která může za velice spe­ci­fic­kých si­tu­ací vést k dvoj­ná­sob­nému zrych­lení je zpětná ite­race jedné smyčky. V běž­ných si­tu­a­cí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ší op­ti­ma­li­zaci, kterou jsem zkou­šel, byly soft­ware ma­naged bu­f­fers zmi­ňo­vané zdezde. Ty fun­gují tak, že data ne­za­pi­suji přímo do paměti, ale do malého bu­f­feru, který se vejde do cache a do paměti ho zko­pí­ruji jen, když se naplní. Tohle při­dalo pár pro­cent výkonu pro řazení intů, ale vedlo to k dras­tic­kému zpo­ma­lení v pří­padě floatů. Vý­sle­dek nebyl nijak pře­svěd­čivý a proto jsem tuhle op­ti­ma­li­zaci na­ko­nec ne­po­u­žil.
  4. To s při­hlád­nu­tím k tomu, že radix sort musí udělat 5 prů­chodů celým polem zna­mená 2.2ns/7taktů CPU na ele­ment v každé ite­raci. Z paměti je přitom třeba číst/za­pi­so­vat 4.8 GB/s (vy­svět­lení 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 ele­ment – pozor jejich čísla nesedí).
  5. V li­te­ra­tuře se uvádí, že radix sort nedělá sek­venční zápis v tom smyslu, že ne­za­pi­suje do jed­noho místa, ale do mnoha, v tomto pří­padě 256. To vede k ne­op­ti­mál­nímu vy­u­žití pro­pust­nosti pamětí a může vést k ne­u­stá­lým vý­pad­kům TLB. Moje ex­pe­ri­menty na mo­der­ním hard­waru však nejsou pře­svěd­čivé aspoň, co se týká TLB. Perf hlásí skoro nulu dTLB-load-missesdTLB-store-misses, ale na druhou stranu cache miss se stále nějaké vy­sky­tují (jako load-miss tak i store-miss).
  6. Se slo­ži­tostí radix sortu je to o něco kom­pli­ko­va­nější. Je pravda, že běží v čase O(kn), ale povaha onoho k je trochu ma­toucí. Když bych měl n růz­ných čísel, nejmenší počet bitů po­třeb­ných k re­pre­zen­taci kaž­dého z nich je log<sub>2</sub>(n). Ono O(kn) je tedy vlastně jen O(n log(n)) v přestro­jení.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články