Budapesti Műszaki és Gazdaságtudományi Egyetem Villamosmérnöki és Informatikai Kar Méréstechnika és Információs Rendszerek Tanszék
Ordering MCMC gyorsítása OpenCL környezetben Önálló laboratórium 1
Készítette Trosztel Mátyás
Konzulens Hajós Gergely
2012. május 29.
Tartalomjegyzék 1. Bevezetés
2
1.1. Az MCMC algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2. Ordering MCMC algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3. Ordering score számítás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.4. Eddigi eredmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.5. OpenCL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2. Tervezés és implementálás
3
2.1. Alapvető technikák . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2. Score számítás gyorsítása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.3. Score indexelés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.4. Score cache mérete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.5. Szülői halmazok generálása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.6. Ordering score számítás gyorsítása . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.7. Validálás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
3. Eredmények
8
3.1. Futási eredmények . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
3.2. További lépések . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
Köszönetnyilvánítás
9
Hivatkozások
10
1
1. 1.1.
Bevezetés Az MCMC algoritmus
Az MCMC (Markov Chain Monte Carlo) algoritmus során egy olyan megfelelően konstruált Markov-láncból mintavételezünk, amely stacionárius eloszlása a céleloszlás. Nagy számú szimulációs lépés esetén a minták a kívánt eloszlást követik. A mintavétel minősége a lépésszám növelésével javítható. Az algoritmus sztochasztikus tulajdonsága miatt több szimuláció eredményét érdemes összevetni a stabilabb adat eléréséért.
1.2.
Ordering MCMC algoritmus
Az Ordering MCMC a sorrendezések tere felett lépked. Ez a tér a struktúrák teréhez hasonlóan szuperexponenciális, de kevésbé „tüskés”, így jobb konvergenciát biztosít. Emellett a sorrendezéssel kompatibilis gráf DAG (irányított körmentes gráf) lesz, ami lehetővé teszi a valószínűségi változók szülői halmazainak egymástól független számítását. Ezt kihasználva az algoritmus párhuzamosítható. A keresési tér csökkentése érdekében maximalizálni szoktuk a szülői változók számát. Ez egy elfogadható közelítés, mert a legtöbb esetben nem áll rendelkezésünk elég adat a nagyobb szülői számmal járó feltételes valószínűségek számításához.
1.3.
Ordering score számítás
Az adatunk valószínűségét egy adott sorrendezésre nézve a következőképp számíthatjuk ki: P (D |≺) =
X
Y
score(Xi , P aG (Xi ) | D),
(1)
G∈Gk,≺ i
ahol score(Xi , P aG (Xi ) | D) annak a valószínűsége, hogy D adat mellett Xi valószínűségi változónak P aG (Xi ) a szülői halmaza. Gk a ≺ sorrendezéssel kompatibilis maximum k darab szülőt tartalmazó gráfok halmaza. P aG (Xi ) az Xi valószínűségi változó szülei a G gráfban. Mivel a szülői halmazok egymástól függetlenül számíthatók: P (D |≺) =
Y X
score(Xi , U | D),
(2)
i U∈Ui,≺
ahol Ui,≺ az Xi valószínűségi változó ≺ sorrendezéssel kompatibilis szülői halmazai. A sorrendezések tere n! ezért zárt formában csak kevés változóra tudnánk kiszámolni, ezért alkalmazzuk az MCMC algoritmust. A Markov-láncot úgy hozzuk létre, hogy a stacionárius eloszlása P (D| ≺) legyen. Ezek után bármilyen f (≺) függvény várható értékét meg tudjuk határozni a következő módon: T 1X E[f | D] ≈ f (≺t ). (3) T t=1
1.4.
Eddigi eredmények
Az előző félévben az MPI (Message Passing Interface) segítségével egy globális gyorsítótárat hoztunk létre, ami lehetővé tette, hogy a párhuzamosan futó szimulációk megosszák egymással a már kiszámolt score-okat, ezzel csökkentve az egy processzre eső score számítást. 2
Eredmények 139 változó esetében: teszt lokális cache elosztott cache (6 klienssel) elosztott cache (12 klienssel)
futásidő 13 568s 5 938s 5 349s
A globális gyorsítótár segítségével több, mint kétszeres gyorsulást sikerült elérni.
1.5.
OpenCL
Az MPI verzió futási idejének fele score számítással telt, ezt szeretnénk most optimalizálni, felgyorsítani. Az MPI azonos gépen lévő processzek esetén nem hatékony megoldás, ezért OpenCLre váltottunk, részben annak reményében, hogy a GPU erősen párhuzamos számításokra készült architektúráját kihasználjuk. Az OpenCL (Open Computing Language) párhuzamos programozás támogatására egy keretrendszert nyújt, amely heterogén platformokon fut (CPU, GPU, FPGA, . . . ). Programnyelve a C99 szabványon alapul, azt a függvényt, ami az OpenCL eszközön fut kernelnek nevezzük. Előnyei: • Az OpenCL programok futásidőben fordulnak, az adott eszköz teljes utasításkészletét ki tudják használni, ezt C-ben csak hardverspecifikusan, többletmunkával tudnánk megtenni. • A futásidőben megadott változókat konstansként tudjuk kezelni, ami további optimalizációt tesz lehetővé, például elágazás és cikluskifejtés. • Vektoros műveletek támogatása. Az optimális kód létrehozását segíti, ezzel hatékonyan ki tudjuk használni az SIMD (single instruction, multiple data) egységeket. • Automatikusan skálázódik a rendelkezésre álló erőforrások szerint. • Szinkronizálás, végrehajtási sorok, események támogatása.
2. 2.1.
Tervezés és implementálás Alapvető technikák
Az egyik alapvető technika a memória elérés optimalizálása. A cache-miss nagyon drága, a gyorsítótárból pár ciklus alatt elérhető az adat, amíg a memóriából csak több, mint 100 ciklus alatt tölthető be. Az összetartozó adatokat tehát tartsuk egymáshoz „közel” a memóriában, és próbáljuk meg a lehető legtöbb műveletet végrehajtani rajta, amíg a cache-ben található. A másik népszerű technika az SIMD (Single Instruction, Multiple Data) egységek kihasználása. Ez SSE (Streaming SIMD Extensions) esetén 8 darab 128 bites regisztert, míg az újabb processzorokon (Intel: Sandy Bridge, AMD: Bulldozer) megtalálható AVX (Advanced Vector Extensions) 16 darab 256 bites regisztert jelent.
3
típus byte short int, float long long, double
SSE 16 8 4 2
AVX 32 16 8 4
A fenti táblázatban látható, hogy milyen típusok esetén egyszerre hány adaton tudjuk elvégezni ugyanezt a műveletet. Az SIMD egységek kihasználásával tehát jelentős gyorsulást érhetünk el. A vektoros utasításkészlet alkalmazásakor nagyon fontos az adatok folytonos tárolása és igazított elérése. Ez azt jelenti, hogy az adatunk memóriacíme az igazítás egész számú többszörösére esik. Vektoros regiszterekbe 16 byte-os igazítással kétszer-háromszor gyorsabban lehet betölteni az adatokat. Gondoskodni kell tehát a megfelelő igazításról, erre a következő segédfüggvényeket használtam: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
void ∗ a l i g n e d _ m a l l o c ( s i z e _ t s i z e , s i z e _ t a l i g n m e n t ) { void ∗ ptr , ∗ a l i g n e d P t r ; p t r = m a l l o c ( s i z e + alignment −1 + s i z e o f ( void ∗ ) ) ; i f ( ! p t r ) return NULL; a l i g n e d P t r = ( void ∗ ) ( ( ( i n t p t r _ t ) p t r + s i z e o f ( void ∗ ) + alignment −1) & ~( al i gnment −1) ) ; ∗ ( ( void ∗ ∗ ) a l i g n e d P t r −1) = p t r ; return a l i g n e d P t r ; } void a l i g n e d _ f r e e ( void ∗ a l i g n e d P t r ) { f r e e ( ∗ ( ( void ∗ ∗ ) a l i g n e d P t r −1) ) ; }
Annyi memóriát kérünk, hogy biztosan el tudjuk tolni egy igazított címre, illetve egy mutatónyi területet fenn kell hogy tartsunk adminisztrációs célra, mert csak ezzel a címmel tudjuk felszabadítani a memóriát. GPU esetén 32, 64, 128 byte-os igazított tranzakcióban történik a memória elérés, így itt különösen fontos, hogy az adatok igazítottan és folytonosan helyezkedjenek el, máskülönben jelentősen lecsökken az elérhető sávszélesség. Véletlenszerű elérés esetén az olvasás 16-szor lesz lassabb, mint optimális esetben. GPU esetén a ciklusok és elágazások is érezhető lassuláshoz vezethetnek. A divergens szálak jelentősen csökkentik a teljesítményt, mert azokat a vezérlőnek sorosítania kell.
2.2.
Score számítás gyorsítása
A score számítás egyszerű műveletek sorozata, amit jól tudunk vektoros műveletekké alakítani. A vektorokban az adott valószínűségi változóhoz tartozó adatok fognak szerepelni, így egy sor áráért 8, illetve AVX támogatás esetén 16 sort tudunk feldolgozni. GPU esetén azt is biztosítjuk ezzel, hogy igazított 32 byte-os tranzakciókban érjük el a memóriát, így a sávszélesség optimálisan lesz kihasználva. Ehhez a bemeneti adatfájlt transzponálni kell.
d1,1 d2,1 d3,1 ...
d1,2 d2,2 d3,2 ...
d1,3 d2,3 d3,3 ...
dnV ariables,1 dnV ariables,2 dnV ariables,3 4
··· ··· ··· ··· ···
d1,nSamples d2,nSamples d3,nSamples ... dnV ariables,nSamples
padding padding padding padding padding
Gondoskodnunk kell a valószínűségi változók első adataihoz tartozó memóriacímek helyes igazításáról (d1,1 , d2,1 , . . . , dnV ariables,1 ). Ehhez igazított allokációt kell használnunk, és a sorokat ki kell egészíteni úgy, hogy a mérete az igazítás egész számú többszöröse legyen. Azt, hogy az adattömbünk mindig elérhető legyen a cache-ben úgy tudjuk biztosítani, hogy a score tábla generálást még az MCMC szimuláció megkezdése előtt elvégezzük. Ha ezt nem tudjuk megoldani, például a hatalmas tár és időigény miatt (lásd 2.4), akkor gyűjtsünk össze annyi hiányzó score-t, amennyit csak tudunk, majd egy lépésben számoljuk ki őket. A feldolgozás sebességén még gyorsíthatunk, ha minden lehetséges szülői halmazszámra külön kernelt fordítunk. Ekkor ugyanis a fordító képes a belső ciklust kifejteni, ami GPU esetén jelentős gyorsulást eredményez.
2.3.
Score indexelés
A már kiszámolt score-ok tárolásához és gyors előhívásához le kell képeznünk egy {változó, szülőhalmaz} konfigurációt egy egyedi indexre, lehetőleg folytonosan, annak érdekében, hogy minél kisebb méretű tömböt tudjunk címezni vele. Szerencsére ezt a feladatot meg lehet oldani optimális módon az l-binomiális számábrázolást felhasználva. Minden egész szám egyértelműen felbontható binomiális együtthatók összegére. Például a 27 felbontása l = 4 esetén: !
!
!
!
6 5 2 1 + + + . 4 3 2 1
27 =
Ezt felhasználva a következőképpen tudjuk egyértelműen leképezni az {Xi , P a(Xi )} konfigurációt: !
n−1 index = Xi ∗ + 4
"
!
!
!
P a(Xi )1 P a(Xi )2 P a(Xi )3 P a(Xi )4 + + + 4 3 2 1
!#
,
(4)
ahol P a(Xi )1 > P a(Xi )2 > P a(Xi )3 > P a(Xi )4 . A fenti módszer akkor képez le optimálisan a [0, n ∗ elemeket P a(Xi )j − 1-el helyettesítjük.
(n−1) ] k
intervallumra, ha a P a(Xi )j > Xi
A különböző szülőhalmaz méreteket eltolással tudjuk egy tömbben tárolni: index = Xi ∗
maxP arents X k=0
!
n−1 + k
|P a(Xi )|−1
X k=0
!
|P a(Xi )|
n−1 + k
X k=0
!
P a(Xi )k , |P a(Xi )| − k
(5)
ahol n a változók száma.
2.4.
Score cache mérete
Az alábbi táblázatban a tábla mérete és a feltöltéséhez szükséges idő található maximálisan három illetve négy szülőszám esetén:
5
változószám 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300
max|P a| = 3 tárigény időigény 5240 B 20.96 ms 90 kB 371.52 ms 479 kB 1.96 s 1550 kB 6.35 s 3838 kB 15.72 s 8034 kB 32.91 s 14 MB 61.39 s 25 MB 105.27 s 40 MB 169.30 s 61 MB 258.88 s 90 MB 6.33 min 128 MB 8.99 min 177 MB 12.41 min 239 MB 16.71 min 315 MB 22.06 min 408 MB 28.59 min 521 MB 36.48 min 656 MB 45.89 min 815 MB 57.02 min 1002 MB 70.06 min 1219 MB 85.22 min 1469 MB 102.71 min 1756 MB 122.77 min 2083 MB 145.63 min 2454 MB 171.55 min 2872 MB 200.78 min 3341 MB 233.60 min 3866 MB 270.28 min 4450 MB 5.19 h 5098 MB 5.94 h
max|P a| = 4 tárigény időigény 10 kB 59.11 ms 393 kB 2.32 s 3262 kB 19.21 s 14 MB 84.80 s 44 MB 266.26 s 112 MB 11.26 min 245 MB 24.67 min 483 MB 48.60 min 878 MB 88.29 min 1497 MB 150.50 min 2425 MB 243.69 min 3763 MB 6.30 h 5637 MB 9.44 h 8192 MB 13.72 h 11 GB 19.43 h 15 GB 26.89 h 21 GB 36.49 h 28 GB 48.66 h 37 GB 63.87 h 48 GB 3.44 nap 61 GB 4.40 nap 77 GB 5.56 nap 97 GB 6.96 nap 120 GB 8.61 nap 148 GB 10.58 nap 180 GB 12.88 nap 217 GB 15.57 nap 261 GB 18.69 nap 311 GB 22.29 nap 369 GB 26.42 nap
Látható, hogy amíg maximálisan három szülőszámmal egzakt keresés esetén is belátható időn belül végzünk, addig négy szülő esetén a tárigény jelentősen nő. Erre megoldás lehet az, ha csak három szülőig számolunk egzakt módon, utána pedig közelítő algoritmusokat alkalmazunk. A négy szülős konfigurációk score-jait pedig csak akkor számoljuk ki, amikor szükség van rájuk.
2.5.
Szülői halmazok generálása
A szülői halmazok generálását szeretnénk párhuzamosan elvégezni. Az lenne az ideális, ha csak egy halmaz kezdő és végső indexét kéne megadni. Ez majd jelentősen egyszerűsíti a párhuzamosítást grid környezetben, mert nem terheljük le a hálózatot előre generált halmazok átküldésével. Az algoritmus a fent ismertetett l-binomiális módszer inverze. A fenti példát megfordítva: rendelkezésünkre áll egy m = 27 index és szeretnénk a hozzá tartozó négy elemű halmazt megtudni. Ehhez először megkeressük azt a legnagyobb n számot, amire még n4 ≤ m, ez lesz a halmaz egyik eleme, ezután ezt kivonjuk m-ből és a többi elemre is elvégezzük ugyanezt. Az algoritmus implementációja: 6
1 int c = nVariables ; 2 f o r ( i n t k = n P a r e n t s ; k > 0 ; k−−) { 3 f o r ( c−−; m_choose ( c , k ) > m; c−−) ; 4 m −= m_choose ( c , k ) ; 5 p a r e n t s [ k −1] = c ; 6 }
Az nk számítás gyorsítására előre kiszámított értékeket használunk. A c meghatározása lineáris kereséssel történik, ez lehetne gyorsabban bonyolultabb kereséssel, de kis n mellett nem nyerünk vele semmit a plusz elágazások miatt.
A gyakorlatban a szülőhalmaz generálása úgy történik, hogy beérkezik az Xi valószínűségi változó, a lehetséges szülőváltozók halmaza H = {Y : Xi ≺ Y }, a kezdő- és végső index. Ekkor a fenti módon generált halmazok elemszáma 0 és |H| − 1 között lesz. Ezzel címezve a lehetséges változók tömbjét, megkapjuk a generált halmazt.
2.6.
Ordering score számítás gyorsítása
Láthattuk hogy a valószínűségi változók szülőhalmazai függetlenül számíthatók, így a score számítást párhuzamosan el tudjuk végezni. Az előző részben azt is láttuk, hogy a halmazok generálását is fel tudjuk bontani, így már nem csak a változók mentén párhuzamosíthatunk. Ez hatékonyabb erőforrás kihasználást tesz lehetővé. Egy fontos megállapítás, hogy a flip MCMC operátor esetén nem kell újraszámolnunk az összes változó score-ját, csak ahol változás történt a szülőhalmazokban. Sőt, ha a cache-elt score kompatibilis az új sorrendezéssel, akkor csak azokat a szülői halmazokat kell legenerálni amelyben az új változó szerepel. Ez jelentős gyorsulást eredményez, főleg akkor, ha a nagy szülőszámú konfigurációkra nem végzünk cachelést.
2.7.
Validálás
Az OpenCL-es kód rossz debuggolhatósága miatt úgy ellenőrzöm a helyes működést, hogy ugyanazt a függvényt megírom C-ben és összehasonlítom a kimenetüket. Viszont különböző optimalizáció és beépített függvények (például lgamma()) eltérő pontossága miatt a lebegőpontos számok nem lesznek pontosan egyenlők. Első gondolat a relatív hiba használata: result1 − result2 < , ahol result2 > result1. result2
Ezzel az a probléma, hogy nullához közeli számoknál rosszul működik. Az IEEE szabvány szerint a lebegőpontos számok bitmintája lexikografikus sorrendben kell hogy kövessék egymást, ezért alkalmazhatjuk a következő trükköt: 1 bool a l m o s t E q f ( f l o a t f 1 , f l o a t f 2 , i n t maxUlps ) { 2 return abs ( ∗ ( i n t ∗ )&f 1 − ∗ ( i n t ∗ )&f 2 ) <= maxUlps ) 3 }
Ahol maxUlps f 1 és f 2 között lévő lebegőpontos számok maximális számát jelöli.
7
3.
Eredmények
3.1.
Futási eredmények
Score számítás: eszköz Native C++ (AMD Phenom II X4 @4.0Ghz) Native C++ SSE, Cache opt. (AMD Phenom II X4 @4.0Ghz) OpenCL (AMD Phenom II X4 @4.0Ghz) OpenCL (AMD Radeon HD4850) OpenCL (AMD Radeon HD7950)
sebesség 3 ms/score 140 µs/score 20 µs/score 110 µs/score 24 µs/score
A CPU-n elért 150-szeres (egy magra vetítve 37,5-szörös) gyorsulás biztató eredmény. Ez önmagában kétszeres javulást hoz az MPI-s verzióhoz képest. Az eredményekből az is kiderül, hogy még a csúcskategóriás Radeon HD7950-es videokártya sem képes hozni egy középkategóriás négymagos CPU szintjét. Ez betudható annak, hogy a score számítás alapvetően kevés és egyszerű műveleteken alapul, így a GPU nem tudja hatékonyan átlapolni a magas globális memória késleltetést. Az ordering részét sajnos még nem volt módom tesztelni, mert nem készült el a kód statisztikai része.
3.2.
További lépések
A célunk, hogy grid környezetben magas változószám mellett hatékonyan tudjuk futtatni a szimulációt. A jelenlegi algoritmust MPI-vel kibővítve nagyfokú párhuzamosítást hozhatunk létre a következő módon: 1. Lefoglalunk X gépet, egy processzel. A gépen történő párhuzamosítás hatékonysági okokból OpenCL-el lesz kivitelezve. 2. Felosztjuk a változókat a gépek között redundáns módon (kis cache esetén minden gépen meglesz a teljes gyorsítótár) 3. A redundáns táblákat tartalmazó gépek felosztják egymás között a táblát, kiszámolják a score-okat, majd adatot cserélnek. 4. Minden Markov-lánchoz választunk egy gépet, aki a láncért lesz felelős. 5. Láncért felelős gépek kiszámolják az új sorrendet, meghatározzák melyik változókat kell kiszámolni és a keresési tér bonyolultságával arányosan (ez fontos, ha a négy változós konfigurációkat csak igény esetén számoljuk ki) a szabad gépekre küldik őket. 6. A segéd gépek legenerálják a szülőhalmazokat, az esetleges hiányzó scoreokat kiszámítják majd visszaküldik a score-ok összegét. 7. A láncért felelős gép megfelelő módon feldolgozza a beérkezett score-okat majd jöhet az újabb MCMC lépés. A feladat felosztása hatékonyan elvégezhető, csak indexeket kell átküldenünk a számítást végző gépekre, nem terheljük a hálózatot. A számítási idő is jól meghatározható, így várhatóan nem lesznek hosszabb üresjáratok az adat megérkezésre várva. 8
Köszönetnyilvánítás Köszönettel tartozom konzulensemnek, Hajós Gergelynek, aki a felmerült problémákra mindig készségesen válaszolt, segítségével hozzájárult a tárgy sikeres teljesítéséhez. Külön köszönet illeti Gézsi Andrást, akinek az ordering MCMC implementációját alapul vettem a program elkészítésénél.
9
Hivatkozások [1] Nir Friedman, Daphne Koller, Being Bayesian About Network Structure, 2000. [2] Dirk Husmeier, Learning Bayesian networks with improved MCMC schemes, Biomathematics & Statistics Scotland. [3] Ian Anderson, Combinatorics of Finite Sets, 1987. [4] Agner Fog, Optimizing software in C++, Copenhagen University College of Engineering, 2012. [5] Intel, Writing Optimal OpenCL, R OpenCL SDK version 1.1 Beta, Intel 2011. [6] NVIDIA, OpenCL Best Practices Guide, 2010.
10
