Někdy je nejchytřejší nedělat nic chytrého (další kapitola nekonečného příběhu o optimalizaci)

Nedávno se na Twitteru objevil článek Fast Nearest Neighbor Queries in Haskell, který představuje novou knihovnu, která dokáže za pomoci velice chytrých algoritmů dělat nearest neighbor dotazy rychleji než jakýkoli jiný nástroj na planetě.

Nearest neighbor dotaz je takový, který má najít nejbližší bod z dané množiny bodů v určitém metrickém prostoru. Naivně se dá spočítat tak, že jeden po druhém projde všechny body s tím, že si pamatuje ten nejbližší a jeho vzdálenost1. Chytré řešení všechny body zaindexuje do stromové datové struktury, která nějakým rafinovaným způsobem rekurzivně rozdělí prostor do oblastí, a během hledání většinu těchto regionů přeskočí, protože je garantováno, že v nich obsažené body nemohou být blíže než současný kandidát. Zmíněná knihovna k tomuto účelu používá cover tree.

Článek rozhodně stojí za přečtení, už jen proto, že popisuje, jak psát rychlý numerický kód v Haskellu. Zcela nepřekvapivě jde hlavně o jednu věc: efektivní využití CPU cache ať už pomocí inlinování datových struktur, aby nebylo třeba nahánět zdivočelé pointery, nebo cache oblivious datových struktur jako třeba Van Emde Boas. Výsledky těchto snah jsou shrnuty v několika grafech s výsledky benchmarků.

A právě tady do příběhu vstupuji já.

Když jsem poprvé zíral na výsledná čísla, něco mi v nich nehrálo. Něco mi na nich nesedělo. Jeden benchmark pracoval se 100000 body o 384 dimenzích, ke každému hledal jednoho nejbližšího souseda měřeného euklidovskou vzdáleností a běžel 13.6 minuty. To by na mém stroji odpovídalo víc jak 2600 miliardám taktů procesoru. Naivní řešení by muselo projít všech 100000 bodů a každý porovnat se 100000 ostatními body a při každém porovnání pracovat s 384 prvky vektoru. Dohromady by bylo třeba 3840 miliard operací.

Skoro to vypadalo, že naivní přístup by nemusel být o moc pomalejší.

Bruteforcing for fun and profit

Sedl jsem si tedy ke stroji a za chvíli vyšvihl jednoduchý program v Céčku, který hledá nejbližší sousedy nejnaivnější možnou metodou.

double euclidean_distance_sq(double *vecs, int veclen, int a, int b) {
  double d = 0;
  for (int i = 0; i < veclen; i++) {
    double t = vecs[veclen * a + i] - vecs[veclen * b + i];
    d += t * t;
  }
  return d;
}

int nn(double *vecs, int veccnt, int veclen, int idx) {
  int nearest = -1;
  double dist = INFINITY;
  for (int i = 0; i < veccnt; i++) {
    double d = euclidean_distance_sq_vec_fl(vecs, veclen, idx, i);
    if (i != idx && d < dist) {
      nearest = i;
      dist = d;
    }
  }
  return nearest;
}

struct idxdist {
  int idx; // -1
  double dist; // INFINITY
};

int veccnt = atoi(argv[1]); // počet vektorů
int veclen = atoi(argv[2]); // délka jednoho vektoru

double *vecs = loadVectorsOrGenerateThemOrWhatever(veccnt, veclen);
struct idxdist *idxdists = initializeIdxsToMinusOneAndDistToMinusInfinity(veccnt);

for (int i = 0; i < veccnt; i++) {
  idxdists[i].idx = nn(vecs, veccnt, veclen, i);
}

Brnkačka. Program zkompilovaný pomocí gcc -O3 -funroll-loops běžel na 10× menších datech 36 vteřin. Aby byl srovnatelný s chytrým algoritmem, musel by běžet něco mezi 8 a 10 vteřinami.

Bylo na čase přitlačit na pilu, na vektorizovanou pilu. Nepodařilo se mi donutit GCC, aby automaticky vektorizoval smyčky a generoval SIMD instrukce7, a proto jsem se musel uchýlit k drastickým opatřením. Našel jsem seznam intrinsics – speciálních funkcí, které GCC přeloží na odpovídající instrukce9 – a ručně jsem vektorizoval funkci, která počítá euklidovskou vzdálenost. Nová verze provádí operace se čtyřmi doubly najednou2,8 a GCC tuto smyčku ještě odroluje, což přidá dalších pár procent výkonu6.

double euclidean_distance_sq_vec(double *vecs, int veclen, int a, int b) {
  __m256d d = _mm256_set1_pd(0.0);
  for (int i = 0; i < veclen; i+=4) {
    __m256d aa = _mm256_load_pd(vecs + veclen * a + i);
    __m256d bb = _mm256_load_pd(vecs + veclen * b + i);
    __m256d t = _mm256_sub_pd(aa, bb);
    d = _mm256_add_pd(d, _mm256_mul_pd(t, t));
  }

  __m256d sum = _mm256_hadd_pd(d, d);
  return ((double*)&sum)[0] + ((double*)&sum)[2];
}

Výsledek není na pohled nijak pěkný, ale mě nešlo o krásu, ale jen a pouze o surovou rychlost. A ta se zvýšila: Vektorizované verzi na výpočet stačí jen 20 vteřin. To je lepší, ale zdaleka ne ideální.

Jako další krok jsem zlepšil využití CPU cache. I když jsou vektory uložené v paměti jeden za druhým a procesor maskuje latenci pamětí načítáním dat dopředu (prefetch), limitujícím faktorem se stala propustnost paměti. Kód počítá vzdálenost jednoho vektoru se všemi ostatními. Onen jeden vektor si drží v cache, a všechny ostatní streamuje z DRAM. Každého se dotkne jednou a pak se hned přesune na další. Když začne porovnávat všechno s dalším vektorem, první vektory jsou už dávno vyhozené z cache a je třeba je opět streamovat z DRAM. To vede k mizernému využití cache.

Řešením je rozdělit pole vektorů do malých bloků, které se vejdou do cache, a počítat podobnost vždycky v jednom z těchto bloků. Když je velikost bloku 64, můžu porovnat 64*64 kombinací, ale z paměti stačí načíst jen 64+64 vektorů.3

int tile = 64;

for (int ti = 0; ti < veccnt; ti += tile) {
  for (int tj = 0; tj < veccnt; tj += tile) {
    int maxi = MIN(ti+tile, veccnt);
    for (int i = ti; i < maxi; i++) { // current vector index
      int maxj = MIN(tj+tile, veccnt);
      for (int j = tj; j < maxj; j++) {

        float d = euclidean_distance_sq_vec_fl(vecs, veclen, i, j);
        if (j != i && d < idxdists[i].dist) {
          idxdists[i].idx = j;
          idxdists[i].dist = d;
        }

      }
    }
  }
}

Kód obsahuje 4 vnořené smyčky, dvě vnější iterují po blocích (tile) a dvě vnitřní iterují obsahem bloků. Tohle uspořádání nevypadá na pohled příliš pěkně, ale plní svůj účel: program běží už jen 8.4 vteřiny a je tedy stejně rychlý jako chytrý algoritmus.

Ale to zdaleka nebylo všechno. V ten okamžik jsem se teprve dostal do ráže a v rukávu jsem měl připraveny ještě dva triky.

Dalším krokem bylo použití float místo double. Float potřebuje 4 bajty namísto osmi, z paměti tedy procesor dokáže protlačit 2× více dat, má 2× lepší využití cache a vektorové instrukce pracují s 8 floaty namísto 4 doublů.

float euclidean_distance_sq_vec_fl(float *vecs, int veclen, int a, int b) {
  __m256 d = _mm256_set1_ps(0.0);
  for (int i = 0; i < veclen; i+=8) {
    __m256 aa = _mm256_load_ps(vecs + veclen * a + i);
    __m256 bb = _mm256_load_ps(vecs + veclen * b + i);
    __m256 t = _mm256_sub_ps(aa, bb);
    d = _mm256_add_ps(d, _mm256_mul_ps(t, t));
  }

  d = _mm256_hadd_ps(d, d);
  d = _mm256_hadd_ps(d, d);
  return ((float*)&d)[0] + ((float*)&d)[4];
}

Tahle funkce je podobná té minulé, jen o něco málo ošklivější a 2× rychlejší. Program s těmito změnami běží jen 4.2 vteřiny – 2× rychleji než chytrý algoritmus.

Všechny tyto změny a úpravy přinesly skoro desetinásobnou akceleraci programu, který sám o sobě nebyl nijak strašně neefektivní. Ale samy o sobě jen optimalizovaly jedno-vláknový program a nesnažily se o nějakou razantní změnu. Poslední logickou metou byla paralelizace. Jak se ukázalo, tato změna byla nejjednodušší a nejrychlejší – stačilo před vnější smyčku přidat

#pragma omp parallel for

a zkompilovat s přepínačem -fopenmp4. Ze smyčky se rázem stala paralelní smyčka, která dokáže do posledního vytížit všechna jádra v systému. Program běžel 1.1 vteřinu. To je podle mě docela dobrý výsledek.

Pro zajímavost paralelní verze běží na mém netbooku s Atomem (2 jádra + hyperthreading, 1.66GHz) 17 vteřin, ±1000 vteřin na první generaci RaspberryPi (1 jádro, 700Mhz) a 51 vteřin na druhé generaci RasberryPi (4 jádra, 900MHz, NEON SIMD). Atom bohužel nepodporuje AVX a tak bylo třeba použít SSSE3, které pracuje jen se 128-bitovými vektory.

Shrnuto do jedné tabulky:

základní verze 36s
vektorizace 20s
lepší využití cache 8.4s
float namísto double 4.2s
paralelizace, OpenMP 1.1s

Kompletní program je k dispozici v tomto gistu.

Slovo na závěr

Na závěr musím dodat pár věcí, které zchladí všechny, kteří už-už vyrážejí slavit do ulic.

Autoři testovali onen chytrý algoritmus na Xeonu taktovaném na 2.8 GHz, takže jejich výsledky jsou o něco málo lepší, než se můžou zdát. Z prezentovaných čísel a grafů se dají těžko vyčíst nějaké absolutní hodnoty, ale vypadá to, že do testů zahrnují jak konstrukci cover tree, tak i jeho dotazování. Z jednoho grafu je vidět, že konstrukce nedominuje celkovému času, ale zabere něco jako 20%-30%. Z jiného se pak dá odvodit, že 50–90% všech přístupů do paměti skončilo jako cache-miss a kolem poloviny všech taktů procesor čekal na paměť a neměl co na práci, zatímco můj naivní program každý takt vykonal 2 instrukce a mnoho z nich pracovalo s osmi floaty najednou. Benchmarky jsem (věrný ideálům diletantství) nesnažil spustit na vlastním stroji, takže je možné, že mi mohlo uniknout něco zásadního.

Navíc: I kdyby byl onen cover tree v některých případech pomalejší, můžou existovat situace, kde bude kralovat5 – jako třeba menší počet dimenzí, jiná metrika nebo větší počet datových bodů. K tomu má stále výhodu parametričnosti, kdy můžu snadno změnit datové typy, reprezentaci a funkci vzdálenosti a nepotřebuje hledat několik bodů najednou, aby využil výhody rozdělení do bloků.

Závěr je následující:

Moderní procesory jsou neuvěřitelně rychlé a někdy hloupé řešení s dobrými konstantami může být rychlejší než chytrý algoritmus, který má na papíře lepší asymptoty. Neslavte však předčasně, vždycky nejdřív musíte vědět, čeho je možné dosáhnout.

Dále k tématu:


  1. Pokud mě zajímá n-nejbližších sousedů, pak je budu udržovat v haldě, která mi umožňuje rychle přistoupit k maximu.
  2. Vektor musí mít délku dělitelnou osmi a funkce neřeší, co se stane, když to neplatí. Není problém to dodělat přidáním malé smyčky nebo switche na konci. Dále je třeba mít na paměti, že _mm256_load_pd potřebuje, aby byla paměť zarovnaná na 32 bajtů, jinak segfaultne. Existuje také _mm256_loadu_pd, která načítá nezarovnanou paměť, ale je o něco pomalejší. V tomto případě je třeba místo malloc použít posix_memalign nebo podobnou funkci, která alokuje zarovnanou paměť.
  3. Velikost bloku musí být zvolena tak, aby se data vešla do nějaké úrovně cache. Pokud je tento parametr dimenzován pro L1 cache, je výsledek rychlejší než když je dimenzován pro L2 cache, ale může pojmout méně dat. Nastavení parametru pro sdílenou L3 cache je ještě pomalejší a má navíc ten problém, že L3 cache je většinou sdílená všemi jádry a každé z nich má k dispozici jen odpovídající část kapacity. Hodnota tohoto parametru je závislá na mnoha parametrech a proto ji není dobré nastavovat napevno.
  4. Více informací o OpenMP se dá najít v knize Programming on Parallel Machines
  5. Tady přichází ke slovu curse of dimensionality, která říká, že čím víc rozměrů má daný prostor, tím jsou všechny úhly bližší pravému úhlu a všechny vzdálenosti jsou si podobnější. To má za následek, že algoritmy jako cover tree a podobné fungují hůře a hůře, protože nemají jak prostor rozdělit na výrazně odlišné oblasti.
  6. I když to nemusí být tak úplně pravda na procesorech, které mají micro-op cache. Na ostatních strojích není nadměrná duplikace kódu také příliš vítaná. Když se program nevejde di instrukční cache, může to mít neblahé následky.
  7. Auto-vektorizace na první pohled vypadá jako triviální optimalizace, kterou kompilátor musí dát na první dobrou. Ale ve skutečnosti to ale není vůbec snadné a často je třeba kompilátoru pomoci. Viz: Auto-vectorization with gcc 4.7.
  8. Paralelismus vektorových instrukcí není to samé, co paralelismus out-of-order superskalárních CPU. OOO paralelismus se pokouší najít různé nezávislé instrukce, které může vykonat najednou, vektorové instrukce naproti tomu provedou stejnou operaci na několika datových položkách najednou. OOO je MIMD, vektorové instrukce jsou SIMD.
  9. Pro všechny, kdo nechtějí používat intrinsics, tu je vector extension GCC.

Flattr this!

This entry was posted in Algo, CPU, Paměť. Bookmark the permalink.

2 Responses to Někdy je nejchytřejší nedělat nic chytrého (další kapitola nekonečného příběhu o optimalizaci)

  1. legoxx says:

    vyborny clanok, dalsie dovody pre keep it simple and stupid je ze inteligentny a prespekulovany kod je casto ovela horsie zrozumitelny a spravovatelny

    • pavel says:

      Nevím proč, ale nemyslím, že tohle byl závěr, který jsi si měl udělat.
      Algoritmizace je universální, autorova optimalizace je hardwarově specifická. Stejně dobře si mohl napsat ad-hoc kód v assembleru.

      Jinak samozřejmě výborný a zajímavý článek. Díky.

Leave a Reply

Your email address will not be published. Required fields are marked *