funkcionálně.cz

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

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

4. 3. 2016 — k47

Ne­dávno se na Twit­teru ob­je­vil článek Fast Nea­rest Ne­i­ghbor Que­ries in Haskell, který před­sta­vuje novou kni­hovnu, která dokáže za pomoci velice chyt­rých al­go­ritmů dělat nea­rest ne­i­ghbor dotazy rych­leji než ja­ký­koli jiný ná­stroj na pla­netě.

Nea­rest ne­i­ghbor dotaz je takový, který má najít nej­bližší bod z dané mno­žiny bodů v ur­či­tém me­t­ric­kém pro­storu. Naivně se dá spo­čí­tat tak, že jeden po druhém projde všechny body s tím, že si pa­ma­tuje ten nej­bližší a jeho vzdá­le­nost1. Chytré řešení všechny body za­in­de­xuje do stro­mové datové struk­tury, která ně­ja­kým ra­fi­no­va­ným způ­so­bem re­kur­zivně roz­dělí pro­stor do ob­lastí, a během hle­dání vět­šinu těchto re­gi­onů pře­skočí, pro­tože je ga­ran­to­váno, že v nich ob­sa­žené body ne­mo­hou být blíže než sou­časný kan­di­dát. Zmí­něná knihovna k tomuto účelu po­u­žívá cover tree.

Článek roz­hodně stojí za pře­čtení, už jen proto, že po­pi­suje, jak psát rychlý nu­me­rický kód v Haskellu. Zcela ne­pře­kva­pivě jde hlavně o jednu věc: efek­tivní vy­u­žití CPU cache ať už pomocí in­li­no­vání da­to­vých struk­tur, aby nebylo třeba na­há­nět zdi­vo­čelé poin­tery, nebo cache ob­li­vi­ous da­to­vých struk­tur jako třeba Van Emde Boas. Vý­sledky těchto snah jsou shr­nuty v ně­ko­lika gra­fech s vý­sledky ben­chmarků.

A právě tady do pří­běhu vstu­puji já.

Když jsem poprvé zíral na vý­sledná čísla, něco mi v nich ne­hrálo. Něco mi na nich ne­se­dělo. Jeden ben­chmark pra­co­val se 100000 body o 384 di­men­zích, ke kaž­dému hledal jed­noho nej­bliž­šího sou­seda mě­ře­ného eu­kli­dov­skou vzdá­le­ností a běžel 13.6 minuty. To by na mém stroji od­po­ví­dalo víc jak 2600 mi­li­ar­dám taktů pro­ce­soru. Naivní řešení by muselo projít všech 100000 bodů a každý po­rov­nat se 100000 ostat­ními body a při každém po­rov­nání pra­co­vat s 384 prvky vek­toru. Do­hro­mady by bylo třeba 3840 mi­li­ard ope­rací.

Skoro to vy­pa­dalo, že naivní pří­stup by ne­mu­sel být o moc po­ma­lejší.

Bru­te­for­cing for fun and profit

Sedl jsem si tedy ke stroji a za chvíli vy­švihl jed­no­du­chý pro­gram v Céčku, který hledá nej­bližší sou­sedy nej­na­iv­nější možnou me­to­dou.

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);
}

Brn­kačka. Pro­gram zkom­pi­lo­vaný pomocí gcc -O3 -funroll-loops běžel na 10x men­ších datech 36 vteřin. Aby byl srov­na­telný s chyt­rým al­go­rit­mem, musel by běžet něco mezi 8 a 10 vte­ři­nami.

Bylo na čase při­tla­čit na pilu, na vek­to­ri­zo­va­nou pilu. Ne­po­da­řilo se mi do­nu­tit GCC, aby au­to­ma­ticky vek­to­ri­zo­val smyčky a ge­ne­ro­val SIMD in­strukce7, a proto jsem se musel uchý­lit k dras­tic­kým opat­ře­ním. Našel jsem seznam in­trin­sics – spe­ci­ál­ních funkcí, které GCC pře­loží na od­po­ví­da­jící in­strukce9 – a ručně jsem vek­to­ri­zo­val funkci, která počítá eu­kli­dov­skou vzdá­le­nost. Nová verze pro­vádí ope­race se čtyřmi doubly na­jed­nou2,8 a GCC tuto smyčku ještě odro­luje, což přidá dal­ších pár pro­cent 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ý­sle­dek není na pohled nijak pěkný, ale mě nešlo o krásu, ale jen a pouze o su­ro­vou rych­lost. A ta se zvý­šila: Vek­to­ri­zo­vané verzi na vý­po­čet stačí jen 20 vteřin. To je lepší, ale zda­leka ne ide­ální.

Jako další krok jsem zlep­šil vy­u­žití CPU cache. I když jsou vek­tory ulo­žené v paměti jeden za druhým a pro­ce­sor mas­kuje la­tenci pamětí na­čí­tá­ním dat do­předu (pre­fetch), li­mi­tu­jí­cím fak­to­rem se stala pro­pust­nost paměti. Kód počítá vzdá­le­nost jed­noho vek­toru se všemi ostat­ními. Onen jeden vektor si drží v cache, a všechny ostatní stre­a­muje z DRAM. Kaž­dého se dotkne jednou a pak se hned pře­sune na další. Když začne po­rov­ná­vat všechno s dalším vek­to­rem, první vek­tory jsou už dávno vy­ho­zené z cache a je třeba je opět stre­a­mo­vat z DRAM. To vede k mi­zer­nému vy­u­žití cache.

Ře­še­ním je roz­dě­lit pole vek­torů do malých bloků, které se vejdou do cache, a po­čí­tat po­dob­nost vždycky v jednom z těchto bloků. Když je ve­li­kost bloku 64, můžu po­rov­nat 64*64 kom­bi­nací, ale z paměti stačí načíst jen 64+64 vek­torů.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 ob­sa­huje 4 vno­řené smyčky, dvě vnější ite­rují po blo­cích (tile) a dvě vnitřní ite­rují ob­sa­hem bloků. Tohle uspo­řá­dání ne­vy­padá na pohled příliš pěkně, ale plní svůj účel: pro­gram běží už jen 8.4 vte­řiny a je tedy stejně rychlý jako chytrý al­go­rit­mus.

Ale to zda­leka nebylo všechno. V ten oka­mžik jsem se teprve dostal do ráže a v rukávu jsem měl při­pra­veny ještě dva triky.

Dalším krokem bylo po­u­žití float místo double. Float po­tře­buje 4 bajty na­místo osmi, z paměti tedy pro­ce­sor dokáže pro­tla­čit 2x více dat, má 2x lepší vy­u­žití cache a vek­to­rové in­strukce pra­cují s 8 floaty na­mí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 po­dobná té minulé, jen o něco málo oš­k­li­vější a 2x rych­lejší. Pro­gram s těmito změ­nami běží jen 4.2 vte­řiny – 2x rych­leji než chytrý al­go­rit­mus.

Všechny tyto změny a úpravy při­nesly skoro de­se­ti­ná­sob­nou ak­ce­le­raci pro­gramu, který sám o sobě nebyl nijak strašně ne­e­fek­tivní. Ale samy o sobě jen op­ti­ma­li­zo­valy jedno-vlák­nový pro­gram a ne­sna­žily se o ně­ja­kou ra­zantní změnu. Po­slední lo­gic­kou metou byla pa­ra­le­li­zace. Jak se uká­zalo, tato změna byla nej­jed­no­dušší a nej­rych­lejší – sta­čilo před vnější smyčku přidat

#pragma omp parallel for

a zkom­pi­lo­vat s pře­pí­na­čem -fopenmp4. Ze smyčky se rázem stala pa­ra­lelní smyčka, která dokáže do po­sled­ního vy­tí­žit všechna jádra v sys­tému. Pro­gram běžel 1.1 vte­řinu. To je podle mě docela dobrý vý­sle­dek.

Pro za­jí­ma­vost pa­ra­lelní verze běží na mém net­booku s Atomem (2 jádra + hy­perthrea­ding, 1.66GHz) 17 vteřin, ±1000 vteřin na první ge­ne­raci RaspberryPi (1 jádro, 700Mhz) a 51 vteřin na druhé ge­ne­raci RasberryPi (4 jádra, 900MHz, NEON SIMD). Atom bo­hu­žel ne­pod­po­ruje AVX a tak bylo třeba použít SSSE3, které pra­cuje jen se 128-bi­to­vými vek­tory.

Shr­nuto do jedné ta­bulky:

zá­kladní verze36s
vek­to­ri­zace20s
lepší vy­u­žití cache8.4s
float na­místo double4.2s
pa­ra­le­li­zace, OpenMP1.1s

Kom­pletní pro­gram je k dis­po­zici 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ž vy­rá­žejí slavit do ulic.

Autoři tes­to­vali onen chytrý al­go­rit­mus na Xeonu tak­to­va­ném na 2.8 GHz, takže jejich vý­sledky jsou o něco málo lepší, než se můžou zdát. Z pre­zen­to­va­ných čísel a grafů se dají těžko vyčíst nějaké ab­so­lutní hod­noty, ale vypadá to, že do testů za­hr­nují jak kon­strukci cover tree, tak i jeho do­ta­zo­vání. Z jed­noho grafu je vidět, že kon­strukce ne­do­mi­nuje cel­ko­vému času, ale zabere něco jako 20%-30%. Z jiného se pak dá od­vo­dit, že 50-90% všech pří­stupů do paměti skon­čilo jako cache-miss a kolem po­lo­viny všech taktů pro­ce­sor čekal na paměť a neměl co na práci, za­tímco můj naivní pro­gram každý takt vy­ko­nal 2 in­strukce a mnoho z nich pra­co­valo s osmi floaty na­jed­nou. Ben­chmarky jsem (věrný ide­á­lům di­le­tant­ství) ne­sna­žil spus­tit na vlast­ním stroji, takže je možné, že mi mohlo unik­nout něco zá­sad­ního.

Navíc: I kdyby byl onen cover tree v ně­kte­rých pří­pa­dech po­ma­lejší, můžou exis­to­vat si­tu­ace, kde bude kra­lo­vat5 – jako třeba menší počet di­menzí, jiná me­t­rika nebo větší počet da­to­vých bodů. K tomu má stále výhodu pa­ra­me­t­rič­nosti, kdy můžu snadno změnit datové typy, re­pre­zen­taci a funkci vzdá­le­nosti a ne­po­tře­buje hledat ně­ko­lik bodů na­jed­nou, aby využil výhody roz­dě­lení do bloků.

Závěr je ná­sle­du­jící:

Mo­derní pro­ce­sory jsou ne­u­vě­ři­telně rychlé a někdy hloupé řešení s dob­rými kon­stan­tami může být rych­lejší než chytrý al­go­rit­mus, který má na papíře lepší asymptoty. Ne­slavte však před­časně, vždycky nejdřív musíte vědět, čeho je možné do­sáh­nout.

Dále k tématu:


  1. Pokud mě zajímá n-nej­bliž­ších sou­sedů, pak je budu udr­žo­vat v haldě, která mi umož­ňuje rychle při­stou­pit k maximu.
  2. Vektor musí mít délku dě­li­tel­nou osmi a funkce neřeší, co se stane, když to ne­platí. Není pro­blém to do­dě­lat při­dá­ním malé smyčky nebo switche na konci. Dále je třeba mít na paměti, že _mm256_load_pd po­tře­buje, aby byla paměť za­rov­naná na 32 bajtů, jinak se­g­faultne. Exis­tuje také _mm256_loadu_pd, která načítá ne­za­rov­na­nou paměť, ale je o něco po­ma­lejší. V tomto pří­padě je třeba místo malloc použít posix_memalign nebo po­dob­nou funkci, která alo­kuje za­rov­na­nou paměť.
  3. Ve­li­kost bloku musí být zvo­lena tak, aby se data vešla do nějaké úrovně cache. Pokud je tento pa­ra­metr di­men­zo­ván pro L1 cache, je vý­sle­dek rych­lejší než když je di­men­zo­ván pro L2 cache, ale může po­jmout méně dat. Na­sta­vení pa­ra­me­tru pro sdí­le­nou L3 cache je ještě po­ma­lejší a má navíc ten pro­blém, že L3 cache je vět­ši­nou sdí­lená všemi jádry a každé z nich má k dis­po­zici jen od­po­ví­da­jící část ka­pa­city. Hod­nota tohoto pa­ra­me­tru je zá­vislá na mnoha pa­ra­me­t­rech a proto ji není dobré na­sta­vo­vat na­pevno.
  4. Více in­for­mací o OpenMP se dá najít v knize Pro­gra­m­ming on Pa­ral­lel Ma­chi­nes
  5. Tady při­chází ke slovu curse of di­mensi­o­na­lity, která říká, že čím víc roz­měrů má daný pro­stor, tím jsou všechny úhly bližší pra­vému úhlu a všechny vzdá­le­nosti jsou si po­dob­nější. To má za ná­sle­dek, že al­go­ritmy jako cover treepo­dobné fun­gují hůře a hůře, pro­tože nemají jak pro­stor roz­dě­lit na vý­razně od­lišné ob­lasti.
  6. I když to nemusí být tak úplně pravda na pro­ce­so­rech, které mají micro-op cache. Na ostat­ních stro­jích není nad­měrná du­pli­kace kódu také příliš vítaná. Když se pro­gram ne­ve­jde di in­strukční cache, může to mít ne­blahé ná­sledky.
  7. Auto-vek­to­ri­zace na první pohled vypadá jako tri­vi­ální op­ti­ma­li­zace, kterou kom­pi­lá­tor musí dát na první dobrou. Ale ve sku­teč­nosti to ale není vůbec snadné a často je třeba kom­pi­lá­toru pomoci. Viz: Auto-vec­to­ri­zation with gcc 4.7.
  8. Pa­ra­le­lis­mus vek­to­ro­vých in­strukcí není to samé, co pa­ra­le­lis­mus out-of-order su­per­ska­lár­ních CPU. OOO pa­ra­le­lis­mus se po­kouší najít různé ne­zá­vislé in­strukce, které může vy­ko­nat na­jed­nou, vek­to­rové in­strukce na­proti tomu pro­ve­dou stej­nou ope­raci na ně­ko­lika da­to­vých po­lož­kách na­jed­nou. OOO je MIMD, vek­to­rové in­strukce jsou SIMD.
  9. Pro všechny, kdo ne­chtějí po­u­ží­vat in­trin­sics, tu je vector ex­tension GCC.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články