YOLO tree

Stará internetová moudrost říká „you only live once“. Život je příliš krátký na jistoty a někdy to prostě musíme risknout. Jako třeba během konstrukce binárního vyhledávacího stromu (BST).

BST je super věc, všichni ho známe, milujeme a napíšeme si jeden nebo dva ještě před snídaní. Ale i přesto má určité vlastnosti, které, i když se zdají naprosto samozřejmé, nemusí zrovna být to, co by člověk chtěl. Konstrukce BST je striktně sekvenční a nedá se (přímo1) naivně paralelizovat. To znamená, že pokud chci nejmenší element nějaké sekvence, musím nejdřív postavit celý BST v čase O(n log(n)) a teprve potom z něj vyhmátnout minimum. V tom okamžiku je ale už pozdě, protože jsem udělal všechnu práci nutnou k seřazení sekvence. V ideálním případě bych udělal jen nezbytně nutné množství práce a konstrukci zbytku bych líně odložil na později.

Naproti tomu je quicksort (QS) nebo radix sort (RS) představuje potenciálně líné algoritmy. Když chci n-tý nejmenší prvek, tak si ho v jednoduše quickselectnu a cestou udělám lineární množství práce vstříc plně seřazenému stupu. Stejně tak jsou všechny QS/RS naivně paralelizovatelé. Jako všechny algoritmy, které vstup rekurzivně rozdělí a pak jednotlivé části nezávisle zpracovávají, se u nich dá triviálně dosáhnout logaritmického paralelního zrychlení.

BST se od QS/RS odlišuje tím, že výsledkem není jen seřazená sekvence, ale dynamická datová struktura, do které můžu přidávat nebo z ní odebírat prvky.

Když chci to nejlepší z obou světů, musím na to jít jako Abe a Aaron z Primeru:

They took from their surroundings what was needed and made of it something more.

Reifikací quicksortu můžu dostat BST, který je dynamický a vyvážený a jeho konstrukce je líná a naivně paralelizovatelná.

Jako při každé reifikaci začnu s nehmotným procesem a zhmotním ho do datové struktury. Když v QS vybírám pivot, reifikuji ho jako uzel stromu s tím, že všechny menší prvky přijdou doleva a ty větší doprava. Takhle postupuji rekurzivně k listům a, když si to moje srdce žádá, strom průběžně vyvažuji.

Strom je tvořen dvěma druhy uzlů – nerealizovanými a realizovanými. Ty nerealizované jsou jednoduché seznamy prvků, realizované uzly jsou pak tvořeny pivotem a dvěma podstromy.

Začnu s nerealizovaným seznamem:

[6, 8, 3, 1, 4, 2, 9, 5, 0, 7]

Vyberu pivot, udělám z něj jeden realizovaný uzel, rozdělím seznam na dva nerealizované seznamy a ty strčím jako podstromy. V tomto případě výběr probíhá v Hoarově stylu první dobrá, ale je možné použít i komplikovanější metody jako best-of-3 nebo quasi best-of-9.

                   (6)
                  /   \
[3, 1, 4, 2, 5, 0]     [8, 9, 7]

(Nerealizované seznamy uvádím vždy v hranatých závorkách. Do kulatých píšu pivot realizovaného uz­lu.)

Takhle můžu pokračovat, dokud budu chtít. Můžu třeba realizovat jen cestu k nejmenšímu prvku v průměrném lineárním čase:

            (6)
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

Jak je vidět, velká část stromu vůbec nebyla vyhodnocena.

Stejně tak můžu strom postavit najednou (podstromy můžu potenciálně stavět paralelně) rekurzivní aplikací tohoto pravidla.

Tato QS inspirovaná konstrukce je úzce spjatá s BST. Knuth o tom psal v The Art of Computer Programming: Když jako pivot používám vždycky první element, je výsledný strom stejný jako když stavím binární strom bez vyvažování. Ve výsledku provedu stejná porovnání, jen v jiném pořadí.

Pokud strom nesplňuje moje požadavky na vyváženost, můžu ho průběžně balancovat.

V mém případě může být levý podstrom příliš velký:

                   (6)
                  /   \
[3, 1, 4, 2, 5, 0]     [8, 9, 7]

Rotace si v tomto případě vynutí částečnou realizaci podstromu.

              (6)
             /   \
          (3)     [8, 9, 7]
         /   \
[1, 0, 2]     (4)
                 \
                  [5]

Po rotaci to bude vypadat takhle:

              (4)
             /   \
          (3)     (6)
         /       /   \
[1, 0, 2]     [5]     [8, 9, 7]

Musel jsme realizovat jednu úroveň stromu. Tento proces může pokračovat rekurzivně až k listům v průměrném lineárním čase (a quicksortovské O(n2) v nejhorším případě).

K vyvažování se nehodí techniky z red-black nebo AVL stromů. Ty potřebují označit každý uzel (buď barvou nebo hloubkou) a to nutně vede k realizaci celého podstromu.

Když jsem o tom přemýšlel, nějakou dobu jsem nevěděl, jak tenhle problém vyřešit. Pak jsem ale v paperu Implementing sets efficiently in a functional language narazil na size balanced trees – tedy stromy vyvažované podle velikosti. Velikost podstromu je lokální vlastnost, kterou znám nehledě na to, zdali je podstrom realizovaný nebo ne. Balancování začne jen, když je jeden podstrom k-krát větší než ten druhý.

Vyvažování ale v principu nemusí být nutné, protože když vybírám pivot náhodně, v průměrném případě je výsledný strom vybalancovaný s průměrnou hloubkou 1.386 log(n). Takže místo balancování bych mohl jen randomizovat pořadí vstupu.

Quicksort styl konstrukce vede k lenosti. Teď jak na dynamické vlastnosti? Chci mít možnost do stromu přidávat a odebírat prvky a zachovat při tom jeho lenost a asymptoty.

Vkládání je triviální. Místo toho, aby byl uzel buď realizovaný nebo ne, může být realizovaný se seznam nerealizovaných změn.

Když například do tohohle stromu

            (6)
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

vložím [ 1.47 4.7 47 ], dostanu

            (6)-[1.47, 4.7, 47]
           /   \
        (3)     [8, 9, 7]
       /   \
    (1)     [4, 5]
   /   \
(0)     [2]

Když budu opět chtít získat minimum, všechny nerealizované seznamy, na které narazím cestou, zrealizuji:

            (6)
           /   \
        (3)     [8, 9, 7, 47]
       /   \
    (1)     [4, 5, 4.7]
   /   \
(0)     [2, 1.47]

V tomto případě jsem jen přidával na konec nerealizovaných seznamů.

Mazání můžu vyřešit dvěma způsoby. Buď striktně nebo líně. Striktní varianta najde mazaný prvek ve stromu a hned ho odstraní. Cestou nutně realizuje všechny uzly (to by mělo být v nejhorším případě O(n) práce). Líná varianta do nerealizovaného seznamu jen přidá tombstone. Když se tento objekt eventuálně dostane do listu, odstaní daný element. Tohle řešení je líné, ale zase může vést k hromadění tombstone objektů v nikdy nerealizovaných částech stromu. To se dá vyřešit tak, že každý podstrom nejen sleduje svou velikost, ale také počet náhrobků.

A to by bylo všechno.


Do tohoto bodu jsem se dostal jenom tak, že jsem přemýšlel, kam až se dá výchozí myšlenka dotáhnout. Pak jsem si ale najednou uvědomil, že jsem už něco podobného viděl.

Fraktální stromy, používané v TokuDB, mají velice podobnou strukturu. Druhá verze fraktálních stromů jsou B-stromy, kde je ke každému vnitřnímu uzlu připnutý buffer akumulující změny. Ty jsou eventuálně (např. když se buffer naplní), propagovány o úroveň níž do bufferů v příslušných uzlech nebo přímo do listů s datovými položkami.

Fraktální stromy jsou podobné mým YOLO stromům, i když k této podobě došly z jiného úhlu. TokuDB se snažilo omezit náhodné IO operace a nahradit je sekvenčním zápisem. Když se buffer naplní, jeho čtení a zápis do nižších úrovní je striktně sekvenční. Já jsem se v případě YOLO stromů zajímal jen o líné zpracování a sledoval ho až do extrému. V obou případech to došlo k podobným závěrům.

Pokud jsou tedy YOLO/fraktal tree reifikacemi quicksortu (resp. multi-way quicksortu), musí existovat podobná konstrukce, která reifikuje merge sort. Ta skutečně existuje a co víc – používá se ve světě databází.

Jde o LSMT – Log Structured Merge Trees v kombinaci s tzv. fractional cascasing – reifikovaný merge sort, který je líný a má ty správné asymptoty jako například O(log(n)) hledání.

Shodou náhod na LSMT byla založena první verze fraktálních stromů z TokuDB.


Dále k tématu:


  1. Můžu paralelně postavit několik BST a pak je naivně sloučit s maximálně logaritmickým zrychlením.

Flattr this!

Posted in Algo | Leave a comment

diff a stromy

Minule jsem tu psal, jak dělat diffy a jak se tato znalost dá využít pro kompresi dat. Teď bych chtěl jít o krok dál k diffování libovolných stromů.

Už tehdy jsem se zmínil, že jsem o tom chvíli přemýšlel, ale jako vždycky jsem se nedostal příliš daleko. V hlavě jsem měl jen řešení několika triviálních případů, ale žádné principy. Pak ale moje obsesivní povaha převzala kontrolu a tendence všechno dotáhnout tak daleko, jak jen to bude možné. Jako vždycky jsem potřeboval algoritmus.

Když se dostanu do problémů, měl bych vždycky se vždycky zeptat WWEKD – What would Edward Kmett do? On má řešení mnoha problémů, tak proč ne zrovna toho současného?

Když jsem se rozhlížel po netu, zjistil jsem, že existují nějaké knihovny, jako třeba diffson, které počítají diffy JSON stromů. Ty ale neumí všechno. diffson dělá diffy jen úroveň po úrovni. Zkontroluje kořen jednoho stromu s druhým a když se liší, začne spolu porovnávat vnořené uzly. Neumí uspokojivě spočítat diff pro rozsáhlejší změny struktury stromu.

Například mám strom

        (5)
       /   \
    (2)     [6]
   /   \
[1]     [3]

a chci zjistit jeho diff se stromem, který jen jinak zcela schodný, ale má jen jednu úroveň navíc.

    (4)
   /   \
[0]     (5)
       /   \
    (2)     [6]
   /   \
[1]     [3]

Výsledný diff by měl být maličký – doslova: sem přidej jeden uzel a hotovo.

A právě tady do hry vstupuje Edward Kmett. Vzpomněl jsem si totiž na jeho přednášku, kde mluvil o succinct datových strukturách a wavelet stromech, které mají efektivní operace rank a select. Pak zmínil, že tohle všechno jde použít k práci s JSON dokumenty bez toho, aby je bylo nutné parsovat. Zmínil, že je (za pomoci succint DS) možné vytvořit velice kompaktní index, s jehož pomocí se dají efektivně číst data z JSON stringu bez jeho naparsování do stromové datové struktury. Hlavní ingrediencí tohoto přístupu byly závorky. Ty určovaly, kde daný podstrom začíná a kde končí. V kombinaci s rychlým rank/select pak můžu efektivně najít konec podstromu a celý ho přeskočit.

Závorky jsou hlavní i pro diffování stromů.

Nejdřív musím strom převést na lineární formu, kde závorky udávají začátek a konec každého z podstromů1.

Z dvou příkladů uvedených výše vzniknou následující sekvence tokenů. Číslo je jeden token stejně jako otevírací a uzavírací závorka.

      (5 (2 1 3) 6)
( 4 0 (5 (2 1 3) 6))

Teď když mám dva stringy, můžu spočítat jejich diff přesně tak, jsem to psal minule. Výsledkem bude něco jako:

+(+4+0             -)

A to je všechno. Diff zahrnuje jen nezbytně nutné množství změn: Přidat jeden uzel na začátek s hodnotou 4, levý podstrom obsahuje nulu, pak následuje zbytek stromu beze změny a nakonec je nový strom ukončen. Změny nejsou tak snadno čitelné jako je JSON diff, ale dokážou zahrnout složité transformace a změny struktury.

Ukážu příklad s rotací stromu. Strom před rotací:

        (6)
       /   \
    (2)     [7]
   /   \
[1]     (4)
       /   \
    [3]     [5]

Po rotaci:

      (4)
     /   \
  (2)     (6)
  / \     / \
[1] [3] [5] [7]

Tokenizace, LCS (nejdelší společná subsekvence) a diff:

A:      (  6 (2 1 ( 4 3       5 ) ) 7)
B:      (4   (2 1     3 ) ( 6 5     7))
LCS:    (    (2 1     3       5     7)
diff:   +4-6     -(-4) +)+(+6  -)-)  +)

Je vidět, že diff obsahuje jen uzly 4 a 6. Všechny ostatní zůstaly na svém místě, protože se nezměnilo jejich relativní pořadí.

Po aplikaci diffu už zbývá jen z nové sekvence tokenů postavit nový strom, který obsahuje aplikované změny.

Tohle všechno se dá dělat líně: Sekvenci tokenů nemusím materializovat, můžu ji produkovat líně, stejně tak i patchování může být líné a výsledek je zase líná sekvence ze které líně stavím strom. Pořád je však třeba postavit zcela nový strom.

Dalším logickým krokem je diff interpretovat jako sérii transformací, kterou je možné provést in-place. Ale to až někdy příště, protože zatím nevím, jak to udělat.


Dodatky:

Běžný algoritmus počítající LCS běží v kvadratickém čase a potřebuje kvadratický prostor. Pomocí metody čtyř rusů se dá dostat na čas O(n2 / log(n)). Tedy tak se to aspoň všude píše, ale nikde jsem nenašel veřejně dostupnou publikaci, která by vysvětlila detaily. Jiný algoritmus potřebuje jen lineární místo, další zas běží v čase O((n+r)log(n)), kde r je editační vzdálenost. Speciální případy (např. neopakující se tokeny) se dají vyřešit rychleji. Dále existují techniky aproximace LCS, které jsou rychlejí než přesný výpočet.

Když jsem začal hledat, jak dělat přímo diff stromů bez převedení na sekvenci, v A survey of potential XHTML diff/merge algorithms jsem narazil na XyDiff. Ten běží v běží v průměrném lineárním čase, produkuje přibližné (ale dostatečně dobré) výsledky a je založen na hashování. Každý podstrom je zhashován a to umožňuje rychle hledat odpovídající podstromy v obou vstupech. Tohle by se dalo použít pro snížení počtu tokenů i v mém přístupu.

Kód tady.


  1. Když jde o strom, který má fixní počet potomků, uzavírací závorky nejsou potřeba, protože jsou implicitní.

Flattr this!

Posted in Algo | Leave a comment

diff a komprese

Někdy se zkrátka daří.

Celé dny se nestane nic a pak přijde den, kdy zjistím, jak funguje unixová utilita diff a použiji těchto znalostí, abych komprimoval 3× rychleji a 70× lépe než ctihodný gzip.

Ale začnu popořadě.

Mám dataset, který je tvořen JSON soubory. Nezáleží kolik přesně jich je, jde o něco jako O(hodně) a v katalogu Dellu jsem nenašel stroj, který by měl RAM tak velikou, že by se do ní všechno vešlo. Nejde jen o náhodné JSON soubory, ale log JSONů (stylem jeden JSON na řádek), které jsou si velice podobné. Záznam může mít něco mezi jedním a sto kilobajty s tím, že se od toho předchozího liší jen změnami na několika málo místech: Nějaké vnitřní vlastnosti mohou být aktualizované, do pole může být přidán jeden nebo dva dílčí objekty a jiný může chybět. Podobnost je veliká a většinou přasahuje 99%.

Když je entropie celého datasetu takto malá, kompresní algoritmy by měly zářit.

Bohužel tomu tak není a Gzip v nejlepším případě komprimoval na čtvrtinu. To mi dlouhou dobu přišlo jako neuspokojivý výsledek s přihlédnutím k obrovské repetici v datech. Dlouho jsem to nijak neřešil, protože disky jsou velké a nestojí skoro nic. Ale data narůstala a narůstala, až do té míry, že jsem o zase začal přemýšlet o lepší kompresi.

Napadlo mě, že když jde o stringy, mohl bych uchovávat jen diffy mezi řádky v logu. To na první pohled vypadá jako jednoduchý problém, který se dá sfouknout za půl hodiny. Když jsem o tom ale začal přemýšlet, došlo mi, že to není tak jednoduché. Naivní algoritmy můžou fungovat v mnoha případech, ale dají se velice snadno zmást a pak nenajdou nejmenší možný diff, ale nějakou monstrozitu, která nic moc neušetří.

Potřeboval jsem něco zásadového, potřeboval jsem algoritmus.

Napadala mě parafráze na Kapitána Willarda z Apokalypsy:

Everyone gets everything he wants. I wanted an algorithm. And for my sins, Wikipedia gave me one. Brought it up to me like room service.

Jak jsem velice rychle zjistil, základem unixové utility diff je algoritmus pro nalezení nejdelší společné podsekvence (longest common subsequence/LCS).

Když mám například dvě sekvence/stringy:

A B D E I J K L X
B C D I J L M N

a správně je zarovnám

A B   D E I J K L     X
  B C D   I J   L M N

je vidět, že jejich LCS je

B D I J L

Když mám tuhle sekvenci, můžu snadno spočítat diff. Z prvního řetězce byly odstraněny znaky A E K X a do druhého přidány C M N. Hotovo.

Řešení LCS se dá najít rekurzivně v nepoužitelně pomalém čase nebo přes dynamické programování s pomocí O(n2) extra prostoru. Paměťové nároky situaci komplikují, protože LCS na úrovni znaků dvou stringů, které mají 100kB, by potřebovalo 40GB pomocnou matici. Z toho důvodu utilita diff nebere v potaz jednotlivé znaky, ale celé řádky, kterých je mnohem méně. To dělá problém lépe zvládnutelný.

Já však neměl text, ale jednořádkové JSONy.

Chvíli jsem začal uvažovat, že bych JSONy naparsoval a počítal jejich diffy na úrovni datových struktur. Ale pak jsem, nejspíš pod vlivem ducha Edwarda Kmetta, dostal nápad. Co kdybych string JSONu rozdělil podle čárek a počítal diff na úrovni těchto segmentů. Velice snadno bych tak dostal bloky textu, které mají v JSONu strukturální význam bez toho, aniž bych ho musel parsovat.

Diff mezi těmito dvěma JSONy

{"a":1,"b":2,"c":3}
{"a":1,"b":47,"c":3}

by vypadal

- 7 6       // odstraní 6 znaků počínaje pozicí 7
+ 7 "b":47, // na pozici 7 vloží nový řetězec

To funguje, ale proto, že JSON obsahuje hodně věcí oddělených čárkou a algoritmus běží v podstatě v kvadratickém čase, není to tak rychlé, jak by to jen mohlo být. Pokud JSON obsahuje vnořené objekty, můžu ho nasekat podle uzavírací složené závorky }1. Těch bývá méně a hledání LCS je pak rychlejší.

Ok, teď tedy vím, jak udělat diff JSONů. Ale je to k něčemu dobré?

Místo vysvětlování ukážu malou tabulku:

jsonlog            2926466947 B   originální soubor
jsonlog.gz          772794031 B   komprimováno gzipem
jsonlog.diffs        80895070 B   komprimováno jako série diffů
jsonlog.diffs.gz     12725864 B   diff + gzip

Jak je vidět, tak v tomto případě je samotné převedení na sérii diffů skoro 10× efektivnější než gzip a kombinace s gzipem vede ke 230× kompresi, velice blízko ke skutečné informační hodnotě datasetu.

Kdyby se někdo zajímal, proč si gzip vede tak zle, tak se musím přiznat, že nevím. Moje teorie je, že nevidí, a tedy ani nedokáže využít, podobnosti ve velkém měřítku a komprimuje v malém pracovním okně. Moje doménové znalosti v tomto případě fungují lépe než obecný gzip2. Uvažoval jsem také o xz. Ten komprimuje o něco lépe, ale je strašlivě pomalý.

Samotné diffování je 3× rychlejší než gzip a diffy se na mém stroji (v jednom vlákně) dekomprimují rychlostí 336 MB/s, což je víc než rychlost čtení/zápisu rotujícího disku. Vyplatí se tedy online dekomprimovat data než je číst z disku.

A to podle mě, stojí za tu námahu.


Zdrojový kód je v tomto gistu a tomto repozitáři


Dodatek #1: Ukazuje se, že zstd s maximální úrovní komprese komprimuje stejně dobře jako můj diff+gzip hack a dekomprese je velice rychlá.

Dodatek #2: Po chvíli testování zstd vypadá jako ideální program na obecnou komprimaci. Úrovní komprese je velice blízko xz/lzma, ale pracuje mnohem rychleji a dokáže data dekomprimovat rychlostí přes 1GB za sekundu.


Poznámky:

  1. Pokud na to chcete jít o něco chytřeji a vyhnout se složeným závorkám v řetězcích, stačí udělat dvoustavové parsování: mimo string a ve stringu. Když jsem mimo string, závorky beru v potaz, když jsem ve stringu, přeskakuji je.
  2. Funguje to takhle dobře jen pokud je podobnost mezi sousedícími záznamy veliká. Když by můj dataset byl tvořen zcela náhodnými JSONy, diffová komprese by neušetřila skoro žádné místo, narozdíl do gzipua, který by stále dokázal srazit velokost na něco kolem jedné čtvrtiny.

Flattr this!

Posted in Algo | 8 Comments

Mýtus o O(1) paměti

Vzpomínám si, jak se mě po jedné z přednášek, kde jsem dlouze mluvil o tom, jak drastický má hardware efekt na naše křehké algoritmy a příčetnost, někdo zeptal, jestli to není jen honba po artefaktech a podivnostech současného hardwaru, a tedy v dlouhodobém měřítku zbytečná práce. Odpověděl jsem něco v tom duchu, že uvažování o několikaúrovňové cache a RAM je užitečné, protože všechny současné procesory – počínaje těmi v mobilech až po ty v obrovských serverech – jsou od sebe skoro k nerozeznání. Možná jsem do plénu hodil zkazku o rychlosti světla a limitaci třírozměrným prostorem, ale nebyl jsem se svojí odpovědí spokojený. Tahle otázka mě donutila provést exkurzi do historie architektury procesorů, ale ani poté jsem neměl skutečnou odpověď.

Tu jsem dostal až teď, když jsem přes Martina Thompsona narazil na čtyřdílný seriál blogpostů The Myth of RAM (2, 3, 4). Jejich autor empiricky i teoreticky ukazuje, že čtení a zápis do paměti – ať už to je cache, RAM, SSD nebo HDD – je O(√N) operace, nikoli O(1). N je v tomto případě velikost paměti, ke které je přistupováno v náhodném a nepředvídatelném pořadí. Nejde nutně o celý dataset, ale jen tu nepředvídatel­nou část.

Tento model je (aspoň podle mého skromného názoru) velký skok kupředu, protože popisuje reálné efekty, která má hardware na chování programů. Ty se s narůstající velikostí dat zpomalují rychleji, než by odpovídalo samotným asymptotám právě proto, že v nich někde figuruje zbloudilá odmocnina N. Například průchod spojovým seznamem není O(N), ale O(N√N). Čtení z hashtabulky není O(1) ale O(√N). Binární hledání není O(log(N)) ale O(log(N)√N). Najednou všechny ad-hoc konstrukce vysvětlující efekty cache v tom kterém případě dostávají uniformitu a řád.

Model rozeznává dva případy: nepředvídatelné a předvídatelné čtení z paměti. Když chci přečíst pole objektů, které leží v paměti jeden vedle druhého, platím √N jen jednou, prefetcher se postará, aby každá další položka dorazila ve skutečném O(1) čase. Prefetcher tedy amortizuje počáteční √N cenu v mnoha následujících předvídatelných operacích.

Články však nezmiňují rozdíl mezi závislými a nezávislými operacemi a vůbec neuvažují paralelní přístup k paměti, MLP a out-of-order superskalární procesory, které zvládají dělat víc věcí najednou včetně několika paralelních dotazů do paměti. Procházení spojovým seznamem je příklad závislých operací. Čtení z hashmapy zase těch nezávislých. Když jsou operace nezávislé, hardware jich může provést několik najednou. Na druhou stranu závislost implikuje sekvenční zpracování (o důsledcích víc tady). O těchto záležitostech se dá uvažovat jako o konfliktu mezi paralelismem a √N cenou paměti. OOO hardware spekuluje hlavně proto, aby mohl začít co nejvíc paměťových transakcí najednou a platit několik √N souběžně.

Celá tahle věc mě nadchla, protože do věci vnáší systém a řád vysvětlující, co se skutečně děje. Není to model, který věci zjednodušuje, ale který je zpřesňuje s ohledem na velikost dat od kilobajtů po petabajty.

Teď mě jen napadá, jestli tento √N faktor nemá nějakou spojitost s cache-oblivious algoritmy, v jichž složitostech skoro vždy figurují právě odmocniny N. To si ale nejspíš nechám jako materiál k přemýšlení pro dlouhé zimní noci v baru HG.

Flattr this!

Posted in CPU, Paměť | Leave a comment

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!

Posted in Algo, řazení | 5 Comments