基于CUDA的多信道鋒電位實時分類方法

蔡瑞初，趙坤垚+，黃禮泊，何 炯，陳 瑤,4

(1.廣東工業(yè)大學 計算機學院，廣東 廣州 510006；2.廣東工業(yè)大學 信息工程學院， 廣東 廣州 510006；3.阿里巴巴達摩院 數(shù)據(jù)庫與存儲實驗室，新加坡 068811; 4.新加坡高等數(shù)字科學中心，新加坡 138602)

0 引 言

在神經(jīng)科學研究中，根據(jù)神經(jīng)元的電位信息來研究神經(jīng)元的編碼機制、分析神經(jīng)元的協(xié)同行為是神經(jīng)科學、生物醫(yī)學和計算機科學等交叉領(lǐng)域科學家使用的基本方法[1]。新型高密度電極探針[2]可同時記錄數(shù)百個密集間隔位點的采集波形，每個位點記錄的波形包含由神經(jīng)元的位置和生理特征等組成的特定信息，將電極探針插入到大腦中并記錄電極尖端周圍的局部電位活動，被觀察到的神經(jīng)元的動作電位稱為鋒電位[3]。由于不同的神經(jīng)元產(chǎn)生的鋒電位不同，對采集到的鋒電位進行提取并將其分配給所屬神經(jīng)元的過程，稱為鋒電位分類[3,4]。

當前，鋒電位分類主要分為低維信道鋒電位分類和高維信道鋒電位分類。針對低維信道分類任務(wù)，文獻[5]基于四元數(shù)的波形特征提取和文獻[5]利用GMM算法去自適應(yīng)分類的方法，都達到了很高的分類準確率，但是兩種方法都是針對低維信道數(shù)據(jù)，數(shù)據(jù)的采樣信道為單信道或者四信道，不能推廣到高維鋒電位數(shù)據(jù)的實時分類應(yīng)用中。隨著高維探針技術(shù)的快速發(fā)展，記錄的高維鋒電位數(shù)據(jù)提供了更多可用于分類的信息，文獻[7]提出的Masked GMM算法，針對GMM算法進行改進，使得高維信道數(shù)據(jù)的計算時間大大縮短，同時也達到了很高的分類準確率，但是其分類算法運行時間過長，影響實時分類場景的應(yīng)用；文獻[8]提出的基于空間和形狀的聚類算法和文獻[9]提出的ISO-SPLIT技術(shù)，在高維鋒電位數(shù)據(jù)的分類準確度和算法執(zhí)行效率方面都有很好的表現(xiàn)，且都達到了實時分類，但是兩個算法分別需要在12個和20個CPU內(nèi)核的工作站上運行，對硬件資源要求較高；文獻[10]采用模板匹配算法并在單GPU上進行了加速，達到了實時分類要求，但是模板匹配方法在鋒電位分類應(yīng)用方面的理論目前還不完善[3]，模板的生成方法對分類結(jié)果影響很大。

基于上述情況，本文提出了基于CUDA的特征掩蔽高斯混合模型，一方面，該模型基于高斯混合模型改進得到，具有很好的抗噪性能，在鋒電位分類中達到了很高的分類精度[7]；另一方面，該模型整個算法流程在CUDA架構(gòu)下進行并行實現(xiàn)并優(yōu)化，大大減少了算法運行時間，同時該算法設(shè)計部署在單GPU上，對硬件要求較低。與文獻[7]相比，本文的并行實現(xiàn)方案在32信道鋒電位數(shù)據(jù)集上，確保分類準確率的同時，實現(xiàn)了170倍的加速，實現(xiàn)了鋒電位的實時分類。

1 背景技術(shù)

1.1 特征掩蔽高斯混合模型

特征掩蔽高斯混合模型(Masked GMM)是由Kadir SN，Goodman DFM等提出來的[11]，該模型主要針對高維聚類分析中存在的樣本特征維度高，從而導致的高斯混合模型出現(xiàn)“維度災(zāi)難”和迭代計算效率低的問題，對于樣本中部分不重要的特征通過設(shè)定特征閾值進行掩蔽，生成掩蔽向量，掩蔽的特征假設(shè)為噪聲且其服從高斯分布，之后將噪聲和未被掩蔽的特征根據(jù)掩蔽向量擬合為新的輸入數(shù)據(jù)，在之后的迭代更新過程中，對于掩蔽特征的更新可以根據(jù)噪聲的均值和方差直接得到，從而達到降低計算量的目的，并極大減少了模型中間運算過程的數(shù)據(jù)存儲空間，因此該算法非常適合高維信道鋒電位數(shù)據(jù)等高維稀疏數(shù)據(jù)的聚類。

該算法主要包括下面6個步驟。

(1)掩蔽向量的生成：訓練樣本集X=[x1,x2,…,xn]T， 樣本數(shù)為N， 維數(shù)為D， 其中Xi=(xi,1,xi,2,…,xi,D)T,xn,i表示第n個樣本點xn的第i個特征。對于所有樣本點的每個特征計算其標準差SD， 同時設(shè)置高低閾值的系數(shù)α與β， 用mn,i表示xn,i的掩蔽向量，mn,i表示如下

(1)

通過掩蔽向量，對于所有樣本的每個特征i， 計算該特征的掩蔽向量為0的特征的均值vi和方差σi2。

(2)使用掩蔽向量重構(gòu)原始數(shù)據(jù)集：對于原始數(shù)據(jù)集，需要根據(jù)掩蔽向量將其擬合為噪聲和真實特征的混合特征，由式(1)可知掩蔽向量mn,i為[0,1]區(qū)間的值，所以將mn,i作為第i個特征的原始數(shù)據(jù)分量的權(quán)重，將1-mn,i作為第i個特征的噪聲分量的權(quán)重，可以計算得到重構(gòu)數(shù)據(jù)集Y， 同時也可以得到重構(gòu)數(shù)據(jù)集的特征方差η， 重構(gòu)數(shù)據(jù)集每個特征yn,i和每個特征的方差ηn,i的計算公式如下

yn,i=mn,ixn,i+(1-mn,i)vi

(2)

zn.i=mn,i(xn,i)2+(1-mn,i)(vi2+σi2)

(3)

ηn,i=zn,i-(yn,i)2

(4)

其中，zn,i表示原始數(shù)據(jù)集的平方通過掩蔽向量擬合后的期望，用于計算特征方差。

(3)M步：計算簇的權(quán)重、均值與協(xié)方差：高斯混合模型通常通過期望最大化算法進行參數(shù)求解，對于特征掩蔽的高斯混合模型，由于掩蔽向量的存在，對于簇k更新權(quán)重ωk、 特征均值μk和協(xié)方差∑k的公式如下

(5)

(6)

(7)

其中，Ck表示簇k中樣本點的集合， |Ck| 表示第k個簇中樣本點的個數(shù)，Mk,i表示第k個簇的第i個特征被掩蔽的樣本點的集合，從式(6)、式(7)可以看到，對于μk和∑k的更新，分為掩蔽特征和非掩蔽特征進行更新，掩蔽特征的更新可以通過噪聲均值和方差直接得到。對于鋒電位分類任務(wù)，鋒電位的非掩蔽特征數(shù)量很少，因此大大降低了存儲空間和計算復雜度。

(4)E步：計算每個樣本點的對數(shù)似然估計：通過上述步驟(3)的計算結(jié)果,可以得到第n個樣本點屬于簇k的對數(shù)似然估計πn,k

(8)

(5)評估聚類結(jié)果并重新分配簇：由于特征掩蔽高斯混合模型在迭代過程中自動確定高斯混合中的聚類數(shù)量，因此需要引入懲罰函數(shù)，該函數(shù)通過抑制具有大量參數(shù)的模型來懲罰過擬合。本文采用貝葉斯信息準則(BIC)

BIC=κln(N)-2ln(L)

(9)

其中，κ為模型中的自由參數(shù)，L為被評估模型的最大似然，N為樣本點數(shù)量。對于特征掩蔽高斯混合模型，若一個簇中有γ個未被掩蔽的特征，那么該簇的自由參數(shù)數(shù)量F(γ) 為

(10)

因此對于有K個簇的掩蔽高斯混合模型，自由參數(shù)κ表示如下

(11)

通過式(8)計算出的πn，k， 將每個樣本點重新分配到似然函數(shù)最大的簇中。對于重新分配后的所有簇，其得分S為

(12)

當?shù)梅諷的變化值小于設(shè)定的閾值則結(jié)束迭代。

(6)簇的刪除和分裂：算法在初始時給定一個可能的最大的簇數(shù)量，之后在M步和E步迭代過程中，通過判斷當前簇是否滿足分裂或刪除的條件，動態(tài)調(diào)整簇的數(shù)量。對于當前的每個簇，如果將其中的樣本點重新分配到對數(shù)似然次大的簇中，聚類評估結(jié)果更優(yōu)，則將聚類評估結(jié)果最優(yōu)的簇刪除，其中的樣本點重新分配到對數(shù)似然次大的簇中。同樣的，對于當前的每個簇，通過M步和E步將其分裂為兩個新簇，并評估聚類結(jié)果，如果聚類評估結(jié)果更優(yōu)，則將評估結(jié)果最優(yōu)的那個簇重新分裂為兩個簇。

由于該模型通過EM算法迭代求解，所以模型的主要運行時間花費在(3)(4)(5)(6)這4個步驟上，因此加快每一輪迭代的計算時間對于實時性的實現(xiàn)極為重要。

1.2 CUDA統(tǒng)一設(shè)備計算架構(gòu)

CPU具有復雜的控制邏輯和大容量的緩存，適合進行事務(wù)管理，處理分支繁雜的任務(wù)，而GPU專為計算密集型、高度并行化的計算而設(shè)計[12]。在NVIDIA公司推出的CUDA編程模型下，使用簡單且易實現(xiàn)的擴展C語言就能編寫GPU并行程序，大大提高了GPU的可編程性。CUDA 架構(gòu)模型[13]如圖1所示，CUDA利用向量將線程組織成為3個不同的層次：線程、線程塊以及線程網(wǎng)格。每個線程都有唯一的線程編號；每個線程都有一個容量較小但速度快的私有寄存器；每個線程塊則擁有一個共享存儲器，該存儲器對塊內(nèi)的所有線程均可見；網(wǎng)格內(nèi)的所有線程都可讀寫一塊相同的全局存儲器和一個只供讀取的常數(shù)存儲器與紋理存儲器。CUDA的主要運算庫包括CUBLAS、CUSOLVER和THRUST等，這些庫能解決典型的大規(guī)模并行計算問題，開發(fā)人員在運算庫的基礎(chǔ)上可以快速、方便地建立起自己的計算應(yīng)用。

圖1 CUDA架構(gòu)模型

2 基于CUDA的并行Masked GMM算法

通過第一小節(jié)介紹的特征掩蔽高斯混合模型和CUDA統(tǒng)一設(shè)備計算架構(gòu)，本節(jié)針對特征掩蔽高斯混合模型的并行性及其優(yōu)化進行分析。

2.1 Masked GMM算法的并行性分析

針對GMM算法，文獻[14]已經(jīng)進行了基礎(chǔ)的分析和實驗，但是其并沒有進行算法并行分析和優(yōu)化，下面針對應(yīng)用于高維稀疏腦電波數(shù)據(jù)的Masked GMM算法的每個步驟進行并行性分析，并根據(jù)分析做相應(yīng)的實現(xiàn)優(yōu)化。

(1)對于掩蔽向量的計算,首先對于每個特征的標準差的計算相互獨立，可以并行計算，通過得到的SDi， 更新特征閾值，對于xn,i， 通過式(1)可知其計算也可以并行執(zhí)行，因此對于每個xn,i分配一個線程去計算該特征的掩蔽向量mn,i， 同樣對于每個特征的噪聲的均值和方差的計算也采用線程編號和特征編號對應(yīng)的方式并行執(zhí)行計算。

(2)對于原始數(shù)據(jù)集的重構(gòu)，由式(2)、式(3)、式(4)可以得到，每個被重構(gòu)的特征點只和原始特征點的數(shù)據(jù)信息有關(guān)，而與其它的特征點無關(guān)，特征點可以和線程塊中的線程一一對應(yīng)，使得重構(gòu)計算可以并行執(zhí)行。

(3)在聚類開始之前，我們需要初始化鋒電位的初始簇編號，在鋒電位分類的應(yīng)用中，我們無法提前確定簇的個數(shù)，相比K-means的迭代初始化方法，本文通過兩個鋒電位的掩蔽向量的海明距離度量兩個鋒電位的相似性，相近的鋒電位初始化為一個簇。對于任意鋒電位，其和已知簇中鋒電位的掩蔽向量的海明距離可以并行計算，從而避免K-means迭代初始化耗費大量時間。

(4)在聚類迭代過程中，通過M步中的式(5)、式(6)、式(7)，可知對于每個簇的權(quán)重，均值和協(xié)方差的計算只與當前的簇中的鋒電位有關(guān)，因此簇與簇之間的M步計算是獨立的，可以并行執(zhí)行，同樣的對于E步式(8)中對數(shù)似然的計算，也以簇為更新的單位，同樣可以并行執(zhí)行。

(5)根據(jù)對數(shù)似然值重新分配簇，需要快速得到每個鋒電位的最大對數(shù)似然值，同時快速計算每個簇所有鋒電位最大似然值的和，可以利用CUDA常用的規(guī)約算法進行并行查找最大似然的簇，同時利用規(guī)約算法求得所有鋒電位的最大似然的和，并行更新得分。

(6)對于刪除步驟，由于對當前存在的每個簇都嘗試刪除并計算得分，最后選取得分最低且低于當前不刪除簇的情況下的得分，所以可以并行嘗試刪除每個簇，并計算得分，最后再判斷是否刪除并進行并行更新。同理對于簇的分裂步驟，也并行嘗試分裂每個簇，對于分裂過程中的M步和E步的調(diào)用，在上面(4)中已經(jīng)分析可以并行執(zhí)行。

綜上，Masked GMM算法中的各個步驟都可以并行執(zhí)行，適合在GPU上進行實現(xiàn)。

2.2 并行Masked GMM算法的優(yōu)化分析

2.2.1 GPU中數(shù)據(jù)存儲與訪問優(yōu)化

(1)鋒電位數(shù)據(jù)傳輸進GPU端進行存儲時按照列優(yōu)先的方式，這樣每個線程對每個數(shù)據(jù)點的特征訪問可以達到合并訪存，減少內(nèi)存訪問的次數(shù)，從而減少訪問內(nèi)存的時間；

(2)對于噪聲的均值v、 方差μ以及E步和M步迭代更新過程中的簇的均值、協(xié)方差和似然估計值等，由于其數(shù)據(jù)量較小，可以放進共享內(nèi)存中，加快線程對其訪問速度；

(3)由于Mask GMM算法迭代過程中臨時變量很多，每次迭代進行分配和釋放臨時變量會帶來巨大時間開銷，因此對于臨時變量，分配為全局變量，在程序整個運行過程中，只分配和釋放一次，減少內(nèi)存分配的時間開銷。

2.2.2 算法的執(zhí)行流程優(yōu)化

圖2 特征解耦優(yōu)化

(2)在GPU端的未優(yōu)化實現(xiàn)中，并行執(zhí)行完M步的簇權(quán)重、均值和協(xié)方差等的更新，將結(jié)果儲存于全局存儲器，用于E步所有簇的更新，在GPU端的優(yōu)化實現(xiàn)中，當M步并行更新完一個簇的權(quán)重，均值和協(xié)方差之后，直接對該簇進行E步的對數(shù)似然估計的更新，流程優(yōu)化如圖3所示，這樣一方面可以避免中間結(jié)果的巨大存儲空間的分配，另一方面將對全局內(nèi)存的訪問轉(zhuǎn)化為對共享內(nèi)存或寄存器的訪問，減少訪問內(nèi)存的開銷。

圖3 處理流程優(yōu)化

(3)由于在M步和E步的迭代更新過程中，每個簇的參數(shù)更新是獨立的，因此我們采用CUDA stream技術(shù)，如圖4所示，其中Copy Engine負責數(shù)據(jù)的傳輸控制，Kernel Engine負責數(shù)據(jù)的并行計算，在一個CUDA stream中可以保證運算按既定的順序執(zhí)行，在不同的stream中的運算可以交疊執(zhí)行，這樣一方面增加了運行的并行性，另一方面可以很好隱藏數(shù)據(jù)傳輸。

圖4 CUDA stream 優(yōu)化

2.2.3 默認參數(shù)的優(yōu)化

由于Masked GMM算法的初始簇為一個比較大的值，默認為500，通過不停的迭代去刪除簇，從而達到最終簇的數(shù)目，通過式(6)、式(7)、式(8)可以看到其更新的計算量和當前存在的簇的數(shù)量有很大的關(guān)系，所以可以在啟發(fā)式的初始化過程中將簇的數(shù)量初始化為一個較小的值，該值隨著初始鋒電位的數(shù)目動態(tài)設(shè)定，初始鋒電位數(shù)量越多，該值越大，然后執(zhí)行迭代，這樣一方面可以減小初始存儲空間的分配，另一方面大大加速算法的收斂速度，減少了計算時間。通過實驗，在不影響分類精度的情況下，算法的參數(shù)優(yōu)化達到了10倍左右的加速。

3 實驗結(jié)果和分析

3.1 實驗平臺和實驗數(shù)據(jù)集

為了驗證基于CUDA加速的Masked GMM算法的性能，本文使用兩個32信道的鋒電位數(shù)據(jù)集進行實驗，數(shù)據(jù)集1在文獻[7]中使用，是由已知的睡眠大鼠的前鹵軀體感覺皮層鋒電位數(shù)據(jù)混合而成，數(shù)據(jù)集2是Kampff實驗室的公開數(shù)據(jù)集，是由麻醉大鼠的初級運動皮質(zhì)(M1)和海馬CA1中記錄得到。本文將每個數(shù)據(jù)集切分為10 s，30 s，60 s和120 s的4個子數(shù)據(jù)集分別進行實驗，并與CPU端的串行實現(xiàn)進行對比評估算法準確性和加速比。

實驗參數(shù)設(shè)置為：CUDA stream數(shù)量為3，初始簇的數(shù)量為500，參數(shù)優(yōu)化后的初始簇數(shù)量為50，高閾值系數(shù)β為4.5，低閾值系數(shù)α為2.0，最大迭代次數(shù)為1000，改變所屬簇的鋒電位數(shù)量閾值為0.05，得分改變閾值為0.01，算法迭代結(jié)束的條件為下列情形之一：①算法運行超過最大迭代次數(shù)1000；②算法在該輪改變所屬簇的鋒電位數(shù)量小于總鋒電位數(shù)量*0.05；③算法兩輪迭代得分差值小于0.01。

實驗所采用的計算平臺見表1。

表1 計算平臺

實驗所采用的數(shù)據(jù)集鋒電位數(shù)量見表2。

表2 數(shù)據(jù)集中的鋒電位數(shù)目

3.2 實驗結(jié)果對比

3.2.1 迭代誤差評估

由于GPU和CPU計算精度的不同，因此Masked GMM算法在GPU端和CPU端運行結(jié)果會產(chǎn)生誤差，為了估計不同平臺代碼的運行誤差，這里分別對兩個數(shù)據(jù)集的30 s的子數(shù)據(jù)集的前5次迭代每個簇的權(quán)重、均值和方差的相對誤差均值進行評估，評估代碼為CPU端代碼和未進行參數(shù)優(yōu)化的GPU端代碼，結(jié)果見表3。

由表3可知，算法在30 s子數(shù)據(jù)集上的前5次迭代的權(quán)重、均值和方差的最大相對誤差為0.000 05，算法的平均運行誤差在鋒電位分類任務(wù)中是可以接受的，也說明了算法流程的執(zhí)行準確性。

表3 相對誤差估計

3.2.2 聚類結(jié)果評估

為了評估并行Masked GMM算法的最終聚類結(jié)果的正確性，將GPU端與CPU端的聚類結(jié)果進行比較，同時GPU端代碼根據(jù)初始參數(shù)的不同，會產(chǎn)生兩個對比結(jié)果，其中GPU默認參數(shù)中初始聚類數(shù)目為500，優(yōu)化后的默認參數(shù)中初始聚類數(shù)目為50，其它參數(shù)不變，同時關(guān)于CPU與GPU之間聚類結(jié)果的比較采用了文獻[9]的對比方法，主要對比兩種方法的主要的簇的數(shù)量和簇中波形的匹配。數(shù)據(jù)集1結(jié)果見表4，數(shù)據(jù)集2結(jié)果見表5。

表4 數(shù)據(jù)集1聚類數(shù)量評估

表5 數(shù)據(jù)集2聚類數(shù)量評估

通過表4和表5，可以看到在數(shù)據(jù)集1和數(shù)據(jù)集2上，未進行參數(shù)優(yōu)化的GPU端代碼和CPU端代碼的最終聚類數(shù)目完全一致，同時對于參數(shù)優(yōu)化的GPU代碼，通過匈牙利算法對優(yōu)化后的GPU代碼的每個簇中的鋒電位和CPU代碼的每個簇中的鋒電位進行最大匹配，得到GPU端的多出來的簇是CPU端的某個簇分裂出來的，這在算法中是可以容忍的，因為在人工校驗階段會由人工去分裂或者合并聚類結(jié)果中的相似簇[7-10]。

3.2.3 運行時間評估

為了評估并行Masked GMM算法的加速效果，我們分別在兩個數(shù)據(jù)集的8個子集上運行CPU端和GPU端的代碼，由于參數(shù)優(yōu)化后的GPU端代碼并不影響聚類結(jié)果的準確性，因此這里GPU端為初始簇數(shù)量為50的實驗結(jié)果，數(shù)據(jù)集1的加速比見表6，數(shù)據(jù)集2的加速比見表7。

由表6和表7可知，在數(shù)據(jù)集1中，GPU端的運行時間相對于CPU端，最大加速比達到了173倍，在數(shù)據(jù)集2中，GPU端的運行時間相對于CPU端，最大加速比達到了160倍，同時在兩個數(shù)據(jù)集中，GPU端運行時間均少于數(shù)據(jù)采集時長，達到了鋒電位數(shù)據(jù)的實時分類。

表6 數(shù)據(jù)集1加速比評估

表7 數(shù)據(jù)集2加速比評估

4 結(jié)束語

本文提出了一種基于CUDA的多信道鋒電位實時分類方法。該方法針對高維信道神經(jīng)元鋒電位數(shù)據(jù)高維稀疏的特性，對Masked GMM模型在單GPU設(shè)備上進行并行實現(xiàn)，同時根據(jù)鋒電位數(shù)據(jù)的特征和算法的處理流程進行優(yōu)化，最終在保證分類精度的情況下，實現(xiàn)了鋒電位的實時分類，并降低了實時分類運算所需的設(shè)備要求。實驗結(jié)果表明，在32信道的鋒電位數(shù)據(jù)集上，基于CUDA的多信道鋒電位實時分類方法相比于文獻[7]的CPU串行實現(xiàn)，達到了最高170倍的加速，同時算法的運行達到了實時，極大地方便了科學研究和實際應(yīng)用，同時，本文的算法具有很大的擴展性，隨著鋒電位數(shù)據(jù)量的增加，可以考慮將算法改進為增量式聚類，使得算法可以在當前硬件設(shè)備上運行更大的鋒電位數(shù)據(jù)并實時分類。