CN109460785B - 一种基于gpu加速的锋电位实时分类的方法 - Google Patents

Abstract

本发明提供一种基于GPU加速的锋电位实时分类的方法，通过将Masked EM算法的各个步骤都划分为互不相关的计算任务，从而使每个步骤内部的计算可以并行化并在GPU运行，使GPU的超强并行计算能力得到发挥，同时各个步骤之间的执行调度和迭代控制由CPU来完成，本发明在基于mask EM聚类算法的基础上，利用GPU对算法进行了加速实现，并对算法的各个流程进行了任务的合理划分和线程的合理分配，使得充分利用计算资源的情况下，一方面保证了结果的准确性，另一方面，大大加快了锋电位分类的计算速度，达到了实时的标准，对于神经科学的研究带来了很大的便利。

Description

一种基于GPU加速的锋电位实时分类的方法

技术领域

本发明涉及一种机器学习技术领域，尤其是一种基于GPU加速的锋电位实时分类的方法。

背景技术

锋电位分类是用于分析电生理数据的一类技术，该技术是根据锋电位形状的相似性将锋电位分组为聚类。原则上，每个神经元倾向于激发特定形状的锋电位，所得到的聚类对应于不同的假定神经元的活动，锋电位分类的最终结果是确定哪个锋电位与这些神经元中的哪一个相对应。由于锋电位信号特征的不稳定性和多电极探针技术的快速发展，采集到的锋电位数据有更多特征，同时锋电位分类算法的运行时间也呈现指数增长，难以满足实际的应用需要。

2014年由Shabnam N.Kadir、Dan F.M.Goodman和Kenneth D.Harriss三人针对EM算法进行的改进，提出Masked EM算法，该算法打破了传统高斯混合模型需要预先确定聚类的初始数目的限制，可以在算法运行过程中，根据特定的准则去分裂或者删除现有的簇，这样使得算法在迭代的过程中可以自适应的调节簇的数量，同时对于高维特征的数据进行聚类需要很大的计算量，该算法将数据的特征根据重要程度计算出一个masked向量，根据设定的阈值掩蔽掉大部分的没有明显特征的维度，同时对有特征的维度进行加权，使得高维特征数据的聚类的计算量大大减少，提高聚类效率，但是该算法对于一个32信道的10s的采样数据，运行时间长达670s，同样不能满足实际应用需要。

发明内容

针对现有技术的不足，本发明提供一种基于GPU加速的锋电位实时分类的方法，该方法可以准确的进行锋电位分类，同时使得计算时间达到实时，满足科学研究等的实际应用需要。

本发明的技术方案为：一种基于GPU加速的锋电位实时分类的方法，通过将MaskedEM算法的各个步骤都划分为互不相关的计算任务，从而使每个步骤内部的计算可以并行化并在GPU运行，使GPU的超强并行计算能力得到发挥，同时各个步骤之间的执行调度和迭代控制由CPU来完成。

进一步的，所述的方法中按照每个锋电位的特征进行任务划分，每个锋电位的每个特征被分配一个线程，从而使得Masked EM算法中掩蔽向量(mask向量)的计算，噪声的均值、方差的计算，原始数据的重构的相关任务能够被大规模并行执行，并且相邻线程的计算任务依赖的数据也相邻，达到合并访存，从而充分利用GPU的显存位宽和带宽。

进一步的，所述的方法利用利用共享存储器作为重构后锋电位数据的高速缓存。

进一步的，所述的方法中利用共享存储器作为重构后锋电位数据的高速缓存，是指首先找出并发线程块所依赖的输入数据，然后线程块将其所依赖的区域数据读到线程块的共享存储器中，之后线程块中的线程从共享存储器中获取计算所需数据进行计算。

进一步的，所述方法中将每个簇的均值和协方差计算与其对数似然估计的计算过程组合起来，从而减少全局内存空间使用的同时减少全局访存次数。

进一步的，所述方法中利用共享存储器作为E步迭代后的当前簇的均值，权重和协方差的高速缓存，并立即进行对数似然估计的运算，省去全局内存的I/O，加速迭代。

进一步的，所述方法中由于对于单次迭代中每个簇之间的协方差和对数似然估计的更新的并行性，对于E步和M步的迭代中关于每个簇的更新，采用CUDA的mutil-stream技术，从而使得程序更好的并发执行。

进一步的，所述方法中程序的总体迭代次数控制和E步与M步的更新控制在CPU端执行，其余部分在GPU端执行。

本发明的有益效果为：在基于mask EM聚类算法的基础上，利用GPU对算法进行了加速实现，并对算法的各个流程进行了任务的合理划分和线程的合理分配，使得充分利用计算资源的情况下，一方面保证了结果的准确性，另一方面，大大加快了锋电位分类的计算速度，达到了实时的标准，对于神经科学的研究带来了很大的便利。

附图说明

图1为本发明的流程示意图。

图2为本发明E步和M步分别在CPU端和GPU端执行的流程图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

如图1所示，一种基于GPU加速的锋电位实时分类的方法，通过将Masked EM算法的各个步骤都划分为互不相关的计算任务，从而使每个步骤内部的计算可以并行化并在GPU运行，使GPU的超强并行计算能力得到发挥，同时各个步骤之间的执行调度和迭代控制由CPU来完成，具体如下：

1)、计算每个特征的掩蔽向量，对于训练集X＝[x₁,x₂,...,x_N]^T，其中，N为样本数，T为矩阵或向量转置符号，x_i＝(x_i,1,x_i,2,....,x_i,D)，D为特征维度，x_n,i表示第n个样本点x_n的第i个特征，m_n,i表示x_n,i的掩蔽向量，m_n,i的计算式如下：

其中，sD_i表示数据集X的第i个特征的标准差，α、β分别表示高低阈值的系数；

通过掩蔽向量计算得到每个特征i的噪声的均值v_i和方差

其计算式如下：

其中，

表示特征i被掩蔽的样本数；

对于每个特征的标准差的计算相互独立，对每个特征分配一个线程进行并行计算，通过得到的标准差sD_i，更新特征阈值，同时，对于每个x_n,i分配一个线程去计算该特征的掩蔽向量m_n,i，同样对于每个特征的噪声的均值v_i和方差

的计算也可以并行执行，因此，对每个特征分配一个线程进行并行计算，对于原始数据在GPU上的存储，采用列优先的方式，这样可以达到合并访问存储的效果，减少读取数据的时间；

2)、通过掩蔽向量和原始数据去重新拟合输入数据，对于原始数据集X中的每个特征点x_n,i通过掩蔽向量m_n,i，可得到原始数据集和噪声组合的虚拟分布，虚拟特征期望y_n,i和方差η_n,i，计算式如下：

y_n,i＝m_n,ix_n,i+(1-m_n,i)v_i

η_n,i＝z_n,i-(y_n,i)²；

其中，z_n,i表示第n个数据点的第i个特征的平方的期望，用于计算方差；

由于数据点的拟合是独立的，只与当前数据点的原始数据和当前特征的噪声均值，方差有关，所以对于每一个数据点x_n,i分配一个线程去拟合特征期望y_n,i和方差η_n,i，同时对于所有特征的期望和方差，首先将其读入共享内存中，在接下来的拟合更新中减少全局内存的访问次数，进而减少特征的期望和方差读取的时间，同时将拟合后的数据的期望和方差也进行列优先存储以达到合并访存的效果；

3)、M步和E步：通过汉明距离计算两个样本点掩蔽向量的距离，并将距离相近的样本点初始化为一个簇，计算簇的权重ω_k、均值μ_k与协方差∑_k，更新样本点对于每个簇的似然对数值，在Masked EM模型中，对于簇k更新权重、特征均值和协方差的公式如下：

其中，ω_k表示第k簇的权重，C_k表示第k簇的数据点的集合，N为数据点的总个数，(μ_k)_i表示第k个簇的第i个特征的均值，(μ_k)_j表示第k个簇的第j个特征的均值，M_k,i表示第k个簇中第i个特征被掩蔽的数据点的集合，(Σ_k)_i,j表示第k个簇的第i个特征和第j个特征的协方差，δ_i,j为第i个特征和第j个特征的相关系数；

第n个样本点属于簇k的对数似然估计π_n,k的计算式如下：

其中，D为数据点的特征维度，T为矩阵或向量装置符号，detΣ_k表示第k个簇的协方差矩阵Σ_k的行列式的值，

为对角线校正项，y_n为拟合后的输入数据，μ_k为簇k的均值，

为第k个簇的第i个特征的方差的逆；

初始化过程中，任意两个样本点之间的汉明距离可以并行计算，然后根据距离进行初始化簇分类，对于第k个簇的权重ω_k，均值μ_k，协方差∑_k的计算，只与第k个簇的数据点有关，第K个簇的对数似然估计与其他簇也独立，如图2所示，采用CUDA的mutil-stream技术，通过设置M个stream，每个stream控制K/M个簇的所有更新，达到高层次并行，同时由于掩蔽向量的存在，使得协方差矩阵和逆矩阵的计算变得麻烦，我们将掩蔽特征和非掩蔽特征分开进行计算，对于掩蔽特征更新，即

采用kernel函数实现，每个线程对应一个特征的更新，对于未被掩蔽的特征

采用CUBLAS函数库中的函数实现，两部分分开实现，一方面大大减少了存储空间，另一方面减轻了编程复杂度；

同样对于

矩阵的逆的计算，矩阵行列式的计算和矩阵的乘法计算都根据掩蔽和非掩蔽特征分开计算，而不是更新完所有簇的权重，均值和协方差再更新对数似然估计，这样不仅可以减少内存占用，同时也减少了访存次数，加快了计算速度；

4)、通过惩罚函数评估聚类，采用贝叶斯信息准则(BIC)通过抑制具有大量参数的模型来惩罚过拟合，即BIC＝κln(N)-2ln(L)，式中，κ为模型中自由参数的个数，N为样本点的数量，L为被评估模型的最大似然；

对于Masked EM模型，若一个簇中有γ个未被掩蔽的特征，那么该簇的自由参数数量

因此，对于有K个簇的掩蔽高斯混合模型，当前的模型自由参数κ表示如下：

从式中可以看出，簇K的自由参数κ的计算只与簇中的样本点及其掩蔽向量有关，所以对于每个簇分配一个线程计算，对于自由参数κ的计算，使用规约算法并行求得所有自由参数的和；

5、计算簇的得分并重新分配簇，通过π_n,k将每个样本点重新分配到似然函数最大的簇中，对于重新分配后的所有簇，其得分

对于每个样本点，采用并行规约算法找到其最大对数似然的所属簇，更新所属簇，并通过并行规约算法求得

更新S；

6)、簇的删除和分裂，在初始时给定一个可能的最大的簇数量，之后在M步和E步迭代过程中，通过判断当前簇是否满足分裂或删除的条件，动态的调整簇的数量；

簇的删除判断如下：对每个簇i，将其中的样本点重新分配到次大的似然函数值所对应的簇中，并更新所有簇中的每个样本点的似然函数的值和每个簇的惩罚项的值，更新之后重新计算得分记为S_i，对于每个簇都尝试删除，可以找到最小的S_new＝min(S₁,S₂,....,S_k)，如果S_new<S，则说明删除之后结果更优，所以对得分最小的簇进行删除，将该簇中的样本点重新分配；

簇分裂判断如下：针对每个簇，通过随机初始化为2个簇，并通过M步和E步进行更新，最终将分成的两个簇，一个簇为当前簇号，另一个簇新分配一个簇号，并通过M步和E步更新之后，重新计算所有簇的得分，如果当前得分小于未分裂得分S，则表示当前簇分裂更优，所以将当前簇分裂为两个；

对于簇的删除操作中，由于对当前存在的每个簇都尝试删除并计算得分，最后选取得分最低且低于当前不删除簇的情况下的得分，所以每个簇分配一个线程并行尝试删除每个簇，并计算得分，最后再判断是否删除并进行并行更新；同理对于簇的分裂步骤，也对每个簇分配一个线程并行尝试分裂。

进一步的，所述的方法中按照每个锋电位的特征进行任务划分，每个锋电位的每个特征被分配一个线程，从而使得Masked EM算法中掩蔽向量(mask向量)的计算，噪声的均值、方差的计算，原始数据的重构的相关任务能够被大规模并行执行，并且相邻线程的计算任务依赖的数据也相邻，达到合并访存，从而充分利用GPU的显存位宽和带宽。

进一步的，所述的方法利用利用共享存储器作为重构后锋电位数据的高速缓存。

进一步的，所述的方法中利用共享存储器作为重构后锋电位数据的高速缓存，是指首先找出并发线程块所依赖的输入数据，然后线程块将其所依赖的区域数据读到线程块的共享存储器中，之后线程块中的线程从共享存储器中获取计算所需数据进行计算。

进一步的，所述方法中将每个簇的均值和协方差计算与其对数似然估计的计算过程组合起来，从而减少全局内存空间使用的同时减少全局访存次数。

进一步的，所述方法中利用共享存储器作为E步迭代后的当前簇的均值，权重和协方差的高速缓存，并立即进行对数似然估计的运算，省去全局内存的I/O，加速迭代。

进一步的，所述方法中由于对于单次迭代中每个簇之间的协方差和对数似然估计的更新的并行性，对于E步和M步的迭代中关于每个簇的更新，采用CUDA的mutil-stream技术，从而使得程序更好的并发执行。

进一步的，所述方法中程序的总体迭代次数控制和E步与M步的更新控制在CPU端执行，其余部分在GPU端执行。

上述实施例和说明书中描述的只是说明本发明的原理和最佳实施例，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。

Claims (8)

1.一种基于GPU加速的锋电位实时分类的方法，其特征在于：通过将Masked EM算法的各个步骤都划分为互不相关的计算任务，从而使每个步骤内部的计算可以并行化并在GPU运行，使GPU的超强并行计算能力得到发挥，同时各个步骤之间的执行调度和迭代控制由CPU来完成，具体如下：

1)、计算每个特征的掩蔽向量，对于训练集X＝[x₁,x₂,...,x_N]^T，其中，N为样本数，T为矩阵或向量转置符号，x_i＝(x_i,1,x_i,2,....,x_i,D)，D为维度，x_n,i表示第n个样本点x_n的第i个特征，m_n,i表示x_n,i的掩蔽向量，m_n,i的计算式如下：

其中，sD_i表示数据集X的第i个特征的标准差，α、β分别表示高低阈值的系数；

通过掩蔽向量计算得到每个特征i的噪声的均值v_i和方差

其计算式如下：

其中，

表示特征i被掩蔽的样本数；

对于每个特征的标准差的计算相互独立，对每个特征分配一个线程进行并行计算，通过得到的标准差sD_i，更新特征阈值，同时，对于每个x_n,i分配一个线程去计算该特征的掩蔽向量m_n,i，同样对于每个特征的噪声的均值v_i和方差

的计算也可以并行执行，因此，对每个特征分配一个线程进行并行计算，对于原始数据在GPU上的存储，采用列优先的方式，这样可以达到合并访问存储的效果，减少读取数据的时间；

2)、通过掩蔽向量和原始数据去重新拟合输入数据，对于原始数据集X中的每个特征点x_n,i通过掩蔽向量m_n,i，可得到原始数据集和噪声组合的虚拟分布，虚拟特征期望y_n,i和方差η_n,i，计算公式如下：

y_n,i＝m_n,ix_n,i+(1-m_n,i)v_i

η_n,i＝z_n,i-(y_n,i)²；

其中，z_n,i表示虚拟分布的第n个数据点的第i个特征的平方的期望，用于计算拟合后输入数据的方差；

由于数据点的拟合是独立的，只与当前数据点的原始数据和当前特征的噪声均值，方差有关，所以对于每一个数据点x_n,i分配一个线程去拟合特征期望y_n,i和方差η_n,i，同时对于所有特征的期望和方差，首先将其读入共享内存中，在接下来的拟合更新中减少全局内存的访问次数，进而减少特征的期望和方差读取的时间，同时将拟合后的数据的期望和方差也进行列优先存储以达到合并访存的效果；

3)、M步和E步：通过汉明距离计算两个样本点掩蔽向量的距离，并将距离相近的样本点初始化为一个簇，计算簇的权重ω_k、均值μ_k与协方差∑_k，更新样本点对于每个簇的似然对数值，在Masked EM模型中，对于簇k更新权重、特征均值和协方差的公式如下：

其中，ω_k表示第k簇的权重，C_k表示第k簇的数据点的集合，N为数据点的总个数，(μ_k)_i表示第k个簇的第i个特征的均值，(μ_k)_j表示第k个簇的第j个特征的均值，M_k,i表示第k个簇中第i个特征被掩蔽的数据点的集合，(Σ_k)_i,j表示第k个簇的第i个特征和第j个特征的协方差，δ_i,j为第i个特征和第j个特征的相关系数；

第n个样本点属于簇k的对数似然估计π_n,k的计算式如下：

其中，D为数据点的特征维度，T为矩阵或向量转置符号，detΣ_k表示第k个簇的协方差矩阵Σ_k的行列式的值，

为对角线校正项，y_n为拟合后的输入数据，μ_k为簇k的均值，

为第k个簇的第i个特征的方差的逆；

初始化过程中，任意两个样本点之间的汉明距离可以并行计算，然后根据距离进行初始化簇分类，对于第k个簇的权重ω_k，均值μ_k，协方差∑_k的计算，只与第k个簇的数据点有关，第K个簇的对数似然估计与其它簇也独立，通过将掩蔽特征和非掩蔽特征分开进行计算，对于掩蔽特征更新，即

采用kernel函数实现，每个线程对应一个特征的更新，对于未被掩蔽的特征

采用CUBLAS函数库中的函数实现，两部分分开实现，一方面大大减少了存储空间，另一方面减轻了编程复杂度；

同样对于

矩阵的逆的计算，矩阵行列式的计算和矩阵的乘法计算都根据掩蔽和非掩蔽特征分开计算，而不是更新完所有簇的权重，均值和协方差再更新对数似然估计，这样不仅可以减少内存占用，同时也减少了访存次数，加快了计算速度；

4)、通过惩罚函数评估聚类，采用贝叶斯信息准则(BIC)通过抑制具有大量参数的模型来惩罚过拟合，即BIC＝κln(N)-2ln(L)，式中，κ为模型中自由参数的个数，N为样本点的数量，L为被评估模型的最大似然；

对于Masked EM模型，若一个簇中有γ个未被掩蔽的特征，那么该簇的自由参数数量

因此，对于有K个簇的掩蔽高斯混合模型，当前的模型自由参数κ表示如下：

从式中可以看出，簇K的自由参数κ的计算只与簇中的样本点及其掩蔽向量有关，所以对于每个簇分配一个线程计算，对于自由参数κ的计算，使用规约算法并行求得所有自由参数的和；

5)、计算簇的得分并重新分配簇，通过π_n,k将每个样本点重新分配到似然函数最大的簇中，对于重新分配后的所有簇，其得分

对于每个样本点，采用并行规约算法找到其最大对数似然的所属簇，更新所属簇，并通过并行规约算法求得

更新S；

6)、簇的删除和分裂，在初始时给定一个可能的最大的簇数量，之后在M步和E步迭代过程中，通过判断当前簇是否满足分裂或删除的条件，动态的调整簇的数量；

簇的删除判断如下：对每个簇i，将其中的样本点重新分配到次大的似然函数值所对应的簇中，并更新所有簇中的每个样本点的似然函数的值和每个簇的惩罚项的值，更新之后重新计算得分记为S_i，对于每个簇都尝试删除，可以找到最小的S_new＝min(S₁,S₂,....,S_k)，如果S_new<S，则说明删除之后结果更优，所以对得分最小的簇进行删除，将该簇中的样本点重新分配；

簇分裂判断如下：针对每个簇，通过随机初始化为2个簇，并通过M步和E步进行更新，最终将分成的两个簇，一个簇为当前簇号，另一个簇新分配一个簇号，并通过M步和E步更新之后，重新计算所有簇的得分，如果当前得分小于未分裂得分S，则表示当前簇分裂更优，所以将当前簇分裂为两个；

对于簇的删除操作中，由于对当前存在的每个簇都尝试删除并计算得分，最后选取得分最低且低于当前不删除簇的情况下的得分，所以每个簇分配一个线程并行尝试删除每个簇，并计算得分，最后再判断是否删除并进行并行更新；同理对于簇的分裂步骤，也对每个簇分配一个线程并行尝试分裂。

2.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：各个步骤的划分按照每个锋电位的特征进行任务划分，每个锋电位的每个特征被分配一个线程，从而使得Masked EM算法中掩蔽向量的计算，噪声的均值、方差的计算，原始数据的重构的相关任务能够被大规模并行执行，并且相邻线程的计算任务依赖的数据也相邻，达到合并访存，从而充分利用GPU的显存位宽和带宽。

3.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：利用共享存储器作为重构后锋电位数据的高速缓存。

4.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：通过将每个簇的均值和协方差计算与其对数似然估计的计算过程组合起来，从而减少全局内存空间使用的同时减少全局访存次数。

5.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：利用共享存储器作为E步迭代后的当前簇的均值，权重和协方差的高速缓存，并立即进行对数似然估计的运算，省去全局内存的I/O，加速迭代。

6.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：由于对于单次迭代中每个簇之间的协方差和对数似然估计的更新的并行性，对于E步和M步的迭代中关于每个簇的更新，采用CUDA的mutil-stream技术，从而使得程序更好的并发执行。

7.根据权利要求1所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：程序的总体迭代次数控制和E步与M步的更新控制在CPU端执行，其余部分在GPU端执行。

8.根据权利要求3或5所述的一种基于GPU加速的锋电位实时分类的方法，其特征在于：利用共享存储器作为重构后锋电位数据的高速缓存，是指首先找出并发线程块所依赖的输入数据，然后线程块将其所依赖的区域数据读到线程块的共享存储器中，之后线程块中的线程从共享存储器中获取计算所需数据进行计算。

Images