99.00000000000000009 problémů s floating point čísly

Floating point čísla a jejich IEEE 754 varianta jsou jednou z těch věcí, které mě nikdy nepřestanou fascinovat. Jde o užitečnou věc, která člověka skoro přesvědčí, že všechno bude ok, že svět má svůj řád, že se stačí odevzdat floating vlnám a už se není třeba o nic starat. Hned potom ale narazí na zaokrouhlovací chyby, nekonečna, NaNy, negativní nulu, denormální čísla a hned z toho narkotického blaha vystřízliví.

V tomto duchu začal 13. sraz přátel PHP – historkou o klasickém případě, kdy floaty nepoužívat – k reprezentaci pe­něz.

Floaty pochopitelně nejsou přesné a zaokrouhlovací chyby vznikají i v triviálních případech. Například 0.1 + 0.2 je 0.30000000000000004. Pokud to šlo peníze, najednou jsem z ničeho vytvořil čtyři femtohalíře.

Tyto chyby jsou malé a nenápadné a člověk si je uvědomí až ve chvíli, kdy mu kancelář začnou obléhat auditoři s beranidly a pokřikují něco o tom, že nesedí účetnictví. Kdyby se peníze ukládaly do intů, tak v nejhorším případě částka přeteče do velkého množství záporných peněz a to je chyba, které si nově vzniklý záporný milionář snadno všimne.

Samozřejmě, na peníze a jiné kvantity, které je třeba vést naprosto přesně, se musí použít něco jako BigInteger nebo BigDecimal.

A to je jenom začátek.

Některé jazyky jako třeba PHP nebo JavaScript přidávají svoje speciality, protože v nich celá čísla automaticky přetékají do floatů.

Například tento PHP kód

for ($i = 0; $i != ($i+1); $i++) { echo "lol"; }

vypadá jak nekonečná smyčka, protože se přece nikdy nemůže stát, že $i je stejné jako $i+1. PHP na to má ale jiný názor. Int eventuálně přeteče do double floatu s nižší přesností, které mají krok mezi dvěma sousedními čísly větší než 1, součet nic nezmění a smyčka skončí. V PHP je int jen double, který čeká na svůj okamžik slávy.

Ve slušných jazycích, by int přetekl do záporných hodnot, ale pořád by platilo, že i a i+1 jsou rozdílné.

Pro zajímavost: celá čísla v intervalu [-224, 224] můžou být v single precision floatu uloženy přesně bez zaokrouhlovacích chyb.


Asociativita je užitečná vlastnost, která je stěžejní pro paralelizaci. Tomu si, jak doufám, žádný slušný člověk nedovolí odporovat. Floaty pochopitelně asociativní nejsou.

val ds = Array.fill(10000)(util.Random.nextDouble)
assert(ds.sum - util.Random.shuffle(ds).sum == 0.0)

val ds = Array.fill(10000) {
        val a, b, c = util.Random.nextDouble
        ((a+b)+c) - (a+(b+c))
}
assert(ds.forall(d => d == 0.0))

Obě aserce selžou, protože součty prováděné v různých pořadích, uvedou různé zaokrouhlovací chyby a to vede k různým finálním výsledkům.

Tohle je poměrně důležité, protože to kompilátorům svazuje ruce v optimalizacích, které můžou bezpečně udělat. Změna pořadí operací může vést ke zrychlení, ale taky k odlišnému výsledku.1


IEEE 754 specifikuje bitovou reprezentaci floatů a toho se dá využít k různým trikům. Jedním z nich je možnost transformovat float na int s tím, že vztah mezi floaty je stejný jako mezi nově vytvořenými inty. Tedy, když platí

a < b

tak platí i

floatToInt(a) < floatToInt(b)

a zároveň existuje transformace zpátky

intToFloat(floatToInt(f)) = f

Jde o pár operací, které ve Scale vypadají takto:

def flip(f: Int) = f ^ (-(f >>> 31) & 0x7FFFFFFF)

def floatToInt(f: Float) = flip(Float.floatToRawIntBits(f))
def intToFloat(i: Int)   = Float.intBitsToFloat(flip(i))

Kód je založený na tomto článku, jen jsem ho upravil pro znaménkové inty Javy. V podstatě jde o to, že pokud je float záporný, převrátím všechny bity kromě znaménka. Tím se bitová reprezentace bude chovat stejně jako int.

Když vynesu průběh intů a floatů na intervalu všech bitových hodnot od samých nul nalevo až po samé jedničky napravo, inty mají takovýhle průběh:

Floaty zase takovýto (NaNy by měly být někde za nekonečny, ale ty si teď dovolím ignorovat):

Kód jen překlopí druhou polovinu intervalu.

Tohohle triku se dá využít k jednoduchému řazení floatů radix sortem. Na začátku floaty převedu na inty zachovávající pořadí, rozjedu radix sort pro inty a na konci výsledek přepíšu zpátky na foaty. Existují i přímé způsoby, jak radixsortovat floaty, ale ty jsou komplikovanější.


Dále k tématu


Pozn.

  1. Některé architektury můžou provádět výpočty s vyšší přesností, třeba 80 bitů namísto 64 a výsledek nakonec zaokrouhlit na 64, viz strictpf.

Flattr this!

Posted in low level | 2 Comments

Závislost je špatná (pro vaše programy i pro váš hardware)

Když jsem nedávno civěl do zdrojáků knihovny breeze, narazil jsem na kód pro slakární součin, který vypadal zhruba takhle:

double sum0, sum1, sum2, sum3 = 0.0;

// ...

for(int i = 0; i < length; i += 4) {
        sum0 += a[i+0] * b[i+0];
        sum1 += a[i+1] * b[i+1];
        sum2 += a[i+2] * b[i+2];
        sum3 += a[i+3] * b[i+3];
}

// ...

return sum0 + sum1 + sum2 + sum3;

Jde o tuto smyčku, čtyřikrát odrolovanou (ve zdrojácích je to 8× ale tady uvádím jen 4× kvůli prostoru):

double sum = 0.0;
for (int i = 0; i < length; i++) {
        sum += a[i] * b[i];
}

Zajímavé je, že součet se neprovádí do jedné proměnné, ale paralelně do čtyř s tím, že na konci funkce se částečné součty zkombinují.

Chvíli jsem dumal nad tím, proč je to napsané zrovna takhle a pak mi došlo, že je to vlastně docela chytrá optimalizace, která zkracuje řetězy závislostí mezi jednotlivými operacemi a může tak vést k většímu instrukčnímu paralelismu.

Když vezmu kód, který sčítá do jedné proměnné

double sum = 0.0;

for (int i = 0; i < length; i += 4) {
        sum += a[i+0] * b[i+0];
        sum += a[i+1] * b[i+1];
        sum += a[i+2] * b[i+2];
        sum += a[i+3] * b[i+3];
}

závislosti mezi jednotlivými operacemi vypadají takhle:

Schéma ukazuje dvě iterace odrolované smyčky a je na něm jasně vidět řetěz sčítání, který musí být proveden sekvenčně. Wide issue a out-of-order procesory, které dokážou odbavit několik nezávislých operací v jenom taktu, teoreticky můžou vykonat všechny čtyři násobení a i += 4 paralelně. Součty však musí počítat sekvenčně.

Naproti tomu verze s částečnými součty má tento graf závislostí:

Jak je vidět, součty na sobě vzájemně nezávisí a hardware je může provést paralelně, jakmile je předchozí várka součtů hotová. A co víc, tento kód má jen jeden krok k plné vektorizaci pomocí FMA instrukcí.

Kompilátor tuto optimalizaci ale nemůže provést jen tak. Floating point operace nejsou komutativní nebo asociativní a když změním pořadí v jakém provádím součty, dostanu lehce odlišný výsledek. V tomhle případě to nevadí a explicitně paralelní verze je tedy preferovaná.

To je jen další příklad, že nejlépe se paralelizují nezávislé operace a platí to na všech úrovních – pro paralelizaci na úrovni strojů, vláken nebo instrukcí.


Dále k tématu:

Flattr this!

Posted in CPU | 2 Comments

Od pohledu dobrý, aneb jak najít skoro stejné obrázky mezi dvěma miliony souborů za méně než deset minut

A tenhle znáte: „Proč je Pentium rychlejší než 486? 486 počítá, ale Pentium jen odhaduje.“ Přestože jde o nejapný vtip o chybě v prvních Pentiích z počítačového pravěku, má v sobě zrnko pravdy: Odhad může být mnohem rychlejší než přesný výpočet.

Odhad, pokud chci, aby byl k něčemu dobrý, není vůbec jednoduchou záležitostí. Nestačí jen vyprodukovat nějaké hausnumero, ale přibližnou odpověď, která má jasně definovanou chybu. Mezi algoritmy, které dělají právě tohle patří kromě klasických velkých jmen jako Bloomfilter, HyperLogLog, Count-min Skech, MinHash a SimHash, patří také perceptuální hashe určené k detekci téměř shodných obrázků.

Všechny zmíněné algoritmy se spoléhají na hashovací funkce, ale některé (MinHash, SimHash a perceptuální hashe o kterých to dneska celé je) jinak než je běžně známé.

Klasické kryptografické hashe, jako je MD5, SHA1 a jejich kamarádi, jsou navrženy tak, že změna jednoho bitu vstupu v ideálním případě znamená změnu poloviny bitů výstupu. Nekryptografické hashe, běžně používané v programování, nemají tuto vlastnost. Typicky jde aspoň o univerzální hash nebo o nějaký k-nezávislý hash. Otisk je přesto výsledkem divokých permutací vstupních dat a neměl by zachovávat jejich korelace. Perceptuální hash je přesně opačný – jde o takový otisk obrázku, který zachovává korelace vstupních dat – když jsou si obrázky na vstupu podobné, budou si s největší pravděpodobností podobné i jejich otisky a míra podobnosti obrázků odpovídá podobnosti hashů.

Algoritmus je jednoduchý:

  1. Zmenším obrázek například na 32×32 pixelů (dělám to kvůli zrychlení následujících výpočtů a ne pro redukci vysokých frekvencí)
  2. Převedu na odstíny šedi (opět kvůli zrychlení následujících kro­ků)
  3. Spočítám 32×32 diskrétní kosinovou transformaci (DCT)
  4. Zredukuji DCT a nechám si jen prvních 8×8 hodnot (kromě té úplně první, která by rozhodila průměr). Tyto hodnoty reprezentují nejnižší frekvence obrázku, tedy ty, které popisují hrubou strukturu.
  5. Spočítám průměrnou hodnotu těchto frekvencí
  6. Vytvořím finální hash o délce 64 bitů – daný bit má hodnotu 1, když je korespondující frekvence vyšší než průměr, 0, pokud je menší

Takto získám perceptuální hash obrázku, kterému nevadí změny barevnosti, světlosti (kvůli tomu, že zaznamenávám rozdíl proti průměru) a ani určité úpravy obrázku samotného. Nejnižší frekvence zaznamenávají nejhrubší kontury a rozdělení světel a stínů a jsou proto odolné vůči změnám v detailech.

Když pak chci porovnat, jak jsou si dva obrázky podobné, spočítám Hammingovu vzdálenost mezi jejich perceptuálními hashi. Výsledkem je počet bitů, ve kterých se oba hashe liší. Když se neliší v žádném bitu, jde o dva vizuálně téměř identické obrázky. Možná byl jednou uložen jako jpg a podruhé jako png. Když se liší jen pár bitů (±5) jde nejspíš o jeden a tentýž obrázek s drobnými úpravami, když se liší ve více bitech, jde s největší pravděpodobností o rozdílné obrázky.

Zkusmo jsem si perceptuální hashování napsal v jazyce Go (respektive bezostyšně ukradl přepsal tohle z Javy). Kompletní kód je k vidění v samostatném gistu.

func ImagePHash(img image.Image) uint64 {
  c := make([]float64, 32)
  for i := 1; i < 32; i++ {
    c[i] = 1
  }
  c[0] = 1 / math.Sqrt(2.0)

  // Reduce size and reduce color
  img = grayscale(resize.Resize(32, 32, img, resize.Bilinear))

  vals := make([][]float64, 32)
  for i := range vals {
    vals[i] = make([]float64, 32)
  }

  bounds := img.Bounds()
  for x := 0; x < bounds.Max.X; x++ {
    for y := 0; y < bounds.Max.Y; y++ {
      _, _, b, _ := img.At(x,y).RGBA()
      vals[x][y] = float64(b)
    }
  }

  /* Compute the DCT */
  dctVals := applyDCT(vals, c)

  // Reduce the DCT and compute the average value (excluding the first term)
  total := 0.0

  for x := 0; x < 8; x++ {
    for y := 0; y < 8; y++ {
      total += dctVals[x][y]
    }
  }
  total -= dctVals[0][0]

  avg := total / float64((8 * 8) - 1)

  // Further reduce the DCT and create 64 bit hash
  hash := uint64(0)

  for x := 0; x < 8; x++ {
    for y := 0; y < 8; y++ {
      if !(x == 0 && y == 0) {
        if dctVals[x][y] > avg {
          hash = hash | (1 << (63-uint64(x*8+y)))
        }
      }
    }
  }

  return hash
}

func applyDCT(f [][]float64, c []float64) [][]float64 {
  F := make([][]float64, 32)
  for i := range F {
    F[i] = make([]float64, 32)
  }

  cosines := [32][32]float64 {}
  for a := 0; a < 32 ; a++ {
    for b := 0; b < 32 ; b++ {
      cosines[a][b] = math.Cos((float64(2*a+1) / float64(2*32)) * float64(b) * math.Pi)
    }
  }

  for u := 0; u < 32; u++ {
    for v := 0; v < 32; v++ {
      sum := 0.0
      for i := 0; i < 32; i++ {
        for j := 0; j < 32; j++ {
          sum += cosines[i][u] * cosines[j][v] * f[i][j]
        }
      }
      sum *= (c[u]*c[v] / 4.0)
      F[u][v] = sum
    }
  }
  return F
}

Jak je vidět, jde o pár řádků kódu a výsledky jsou překvapivě přesné.

(nebo, nebo, nebo)

Přesto nejde o neomylný algoritmus. Kód si pohrává s pravděpodobnostmi a někdy, když se plete, bývají výsledky veselé. Zvlášť, když zvýším toleranci podobnosti.

Teď už zbývá jen vyřešit, jak to dělat rychle, když chci projít v titulku zmíněné dva miliony obrázků pod deset minut.

Samotné porovnání p-hashů se může zdát jako dostatečně rychlé. V ideálním případě jde o pouhé 2 instrukce: xor a popcnt + něco málo okolo.1 Přesto, pokud chci porovnat všechny hashe se všemi ostatními, musím udělat kvadratické množství práce a velké číslo se umocněním na druhou bohužel nezmenší.

Naštěstí můžu všechno znatelně urychlit pomocí postupu známého jako locality-sensitive hashing (LSH) – další z metod, která si hraje s pravděpodobnostmi a obětuje něco málo přesnosti ve prospěch masivního zrychlení.

LSH využívá hashe zachovávající lokalitu (podobné věci mají podobné hashe) a snaží se položky rozdělit do skupin, kdy v každé skupině jsou položky, které si můžou být podobné, ale nejsou tam ty, které jsou očividně rozdílné.

V tomto případě zkonstruuji LSH pro Hammingovu vzdálenost tak, že 64 bitů hashe rozsekám na 8 pruhů (band) po 8 bitech a každý hash přiřadím do skupiny (bucket) určenou bity každého pruhu.

To znamená, že když prvních osm bitů hashe má hodnotu 0xFF, v prvním pruhu ho přiřadím do bucketu 0xFF, když má druhých 8 bitů 0x0A, v druhém pruhu spadne do skupiny 0x0A a tak dále.

Celé to funguje na jednoduchém předpokladu: Pokud jsou si dva hashe velice podobné, s největší pravděpodobností se shodují v bitech jednoho nebo více pruhů.2 Když chci najít všechny hashe podobné hashi A, vezmu jeho prvních 8 bitů, podívám se do odpovídajícího bucketu prvního pruhu a spočítám podobnost jen s těmito kandidáty, pak vezmu dalších osm bitů a spočítám podobnost s těmito kandidáty. To udělám jednou pro každý pruh a mám výsledné podobné hashe. Protože jde o pravděpodobnostní metodu, může se stát, že některé skutečně podobné hashe to přeskočí, protože se vždy o jeden bit lišily ve všech pruzích a proto skončily v odlišných bucketech a vůbec nebyly vybráni jako kandidáti. Pokud se zajímám o velice podobné hashe, je tato šance velice malá. Pokud chci v takto vytvořeném LSH najít hashe lišící se v méně než 8 bitech, je garantováno, že všechny dostanu, protože se aspoň v jednom pruhu musí shodovat.

S LSH, kdy mám 8 pruhů po 8 bitech, mi musím udělat jen 8/256 práce – 32× méně, protože se musím podívat do 8 bucketů (jeden v každém pruhu) a spočítat podobnosti ve všemi kandidáty, kterých je 1/256. Pokud to nestačí, můžu tyto dva parametry zvolit podle toho, jak je třeba: když chci víc přesnosti zvolím menší šířku pruhu (12×5 bitů), pokud chci více rychlosti, zvolím širší pruhy (6×10 bitů).

No a s tímto „zlepšovákem“ dokážu bez větších problémů porovnat ony dva miliony obrázků pod 10 minut.

pattern := "/path/to/images/*"
fs, _ := filepath.Glob(pattern)

files  := make([]string, 0, len(fs))
hashes := make([]uint64, 0, len(fs))

for _, fullPath := range fs {

  reader, err := os.Open(fullPath)
  if err != nil { continue }

  img, _, err := image.Decode(reader)
  if err != nil { continue }

  hash := ImagePHash(img)
  files  = append(files, fullPath)
  hashes = append(hashes, hash)

  reader.Close()
}


// create LSH map

const Bands = 8
const BandBits = 8
const BandMask = (1 << BandBits) - 1

bands := [Bands][BandMask+1][]IdxHashPair {}

for idx, hash := range hashes {
  for b := 0 ; b < Bands ; b++ {
    h := hash << uint32(b*BandBits) & BandMask
    bands[b][h] = append(bands[b][h], IdxHashPair { idx, hash })
  }
}


// compute similarities

for idx, hash := range hashes {
  for b := 0 ; b < Bands ; b++ {
    h := hash << uint32(b*BandBits) & BandMask

    for _, candidate := range bands[b][h] {
      if idx <= candidate.Idx { continue }

      diff := differentBits(hash, candidate.Hash)
      if diff < 6  {
        similarity := 1.0 - (float64(diff) / float64(64))
        fmt.Printf("%s %s %d\n", files[idx], files[candidate.Idx], similarity)
      }
    }
  }
}

To je všechno. Záhada vyřešena.


Protože tuto techniku používám na několika místech (například pro detekci podobných článků na devblozích, tady, na téhle maličkosti a na mnoha dalších místech), tak jsem si proto napsal vlastní knihovnu nazvanou prozaicky sketches, která obsahuje implementaci PHashe, skečů, LSH, ALS, FunkSVD nebo DBSCAN clusterování a dalších nesouvisejících vě­cí.


Dále k tématu:


Pozn:

  1. V Go se bohužel k popcnt bez tanečků v assembleru nedá dostat a tak je třeba si napsat utilitu, která dělá to samé, ale potřebuje o něco víc instrukcí.
  2. Pravděpodobnost, že LSH odhalí dva kandidáty je 1 – (1 – pn)b, kde p je podobnost mezi dvěma kandidáty, n šířka pruhu a b je počet pruhů. Pomocí tohoto vztahu můžu najít hodnoty n a p, které budou fungovat pro požadovaný práh podobnosti.

Flattr this!

Posted in Algo, DS | Leave a comment

Sketches – slajdy

Sketches – slajdy

Flattr this!

Posted in Uncategorized | Leave a comment

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