李小林，王 静，张元孜，黄世国

(福建农林大学 计算机与信息学院，福州 350002) E-mail：fjhsg25@126.com

1 引 言

帝王蝶算法(Monarch Butterfly Optimization，MBO)是一种新的群体智能优化算法，通过模拟帝王蝶迁徙的行为实现寻优过程[1].为了避免算法过早收敛，通过结合和声搜索操作[2]或结合自适应交叉算子和贪婪策略[3]对该算法进行了改进.冯艳红等[4]在帝王蝶算法上结合差分变异策略来生成新个体，提高了获得折扣{0，1}KP问题全局最优解的概率.针对TSP问题，Wang等[5]提出了一种二进制帝王蝶优化算法.

定量构效关系(Quantitative Structure-Activity Relationship，QSAR)的目标是建立化合物的结构性质和相应的生物活性之间的关系[6].化合物的结构由各种特征(分子描述符)进行编码来描述.然而，并不是所有的特征都与其生物活性相关.特征太多可能导致模型过拟合，特征太少则可能会导致欠拟合[7].只包含相关和无冗余特征的最优特征子集能够提高QSAR模型预测的精度和可解释性[8].因此，特征选择成为QSAR研究的关键预处理步骤.

元启发式算法因其较强的全局搜索能力和潜力已成功应用于QSAR特征选择.这些代表性的方法包括遗传算法[9，10]，粒子群优化[11，12]，差分进化算法等.

针对QSAR特征选择需在优先满足最大化准确度前提下同时最小化特征数的多目标要求，本研究在MBO算法的粒子排序步骤中采用了非支配排序算法，并修改了调整算子；采用准确度优先的策略，改善了传统多目标算法在低准确度区域耗费大量计算资源的缺陷；采用分小组使用不同突变策略的方法解决了算法早熟的问题.实验结果表明，这些策略的综合应用提升了MBO算法的性能，实现了QSAR建模提高准确度并减少特征数的目的.

2 定量构效建模

为了建立化合物结构特征与其活性之间的关系，采用机器学习方法来建立二者的关联.其中偏最小二乘法(Partial Least Squares，PLS)常用于预测化合物结构特征与其相应的生物活性或化学性质之间的关系.偏最小二乘法通过具有因子分析的线性多元模型将两个数据矩阵x和y联系起来[13].在QSAR建模中，x表示所有化合物的特征向量，y表示相应化合物的活性.

用Q2值衡量QSAR建模的精度，定义如式(1)所示：

(1)

3 帝王蝶优化算法

帝王蝶优化算法模拟了帝王蝶的迁移行为[1].其行为简单描述为：1)帝王蝶的整个种群被分为两个子种群，称为子种群1和子种群2；2)每个子种群中更优的子代个体将分别由迁移操作算子或调整算子产生，并取代其父代个体；3)当前最优的个体自动存活到下一代.

假设问题对应于一个d维解空间，初始化产生一个大小为NP的种群P=(X1，X2，…，XN)，第i个个体在解空间的位置可以表示为一个d维向量Xi=(xi，1，xi，2，…，xi，d)T，显然，每个个体的位置都对应着目标问题的一个潜在解.其迁移算子和调整算子描述如下.

设NP1和NP2分别为两个子种群中个体数量，p为子种群1个体的占比，即NP1=ceil(p×NP)，NP2=NP-NP1.其中，函数ceil(c)表示向上取整.迁移算子的数学模型描述如式(2)所示：

(2)

其中，i=1，2，…，N1，t表示当前迭代次数，xr1和xr2分别表示帝王蝶个体r1和r2的位置，这两只帝王蝶分别随机取自两个子种群.标量r=rand×peri，rand是[0，1]间的一个随机数；peri代表迁移周期.

在调整算子中，若rand≤p，则子种群2中个体按式(3)更新：

(3)

其中，i=1，2，…，N2.xbest表示全局最优个体的位置.

否则，个体将按式(4)更新：

(4)

其中，xr3表示从子种群2中随机选择的一个帝王蝶的位置.

(5)

(6)

其中，α是权重因子，由α=Smax/t2计算所得，Smax表示帝王蝶的最大步长.dx表示帝王蝶的步长，由Levy函数计算.BAR表示调整率.以最小值优化为例，MBO算法的步骤如算法1所示.

算法1.MBO算法伪代码

初始化参数：种群大小NP，子种群1中个体的比例p，调整率BAR，迁移周期peri，权重因子α，最大步长Smax，t=1，最大评价次数T；

随机初始化大小为NP的种群；

计算种群中每个个体的适应度值；

Whilet≤Tdo

根据适应度值对所有个体进行排序；

适应度值较小的前p×NP个个体组成子种群1，其余个体组成子种群2；

For子种群1中的每个个体

根据公式(2)进行个体的更新；

End

For子种群2中的每个个体

根据公式(3-6)进行个体的更新；

End

将两个更新后的子种群合并成一个种群；

计算种群中每个个体的适应度值；

t=t+1；

End

输出最优个体的位置xbest及其对应的适应度值f(xbest).

4 准确度优先的多目标帝王蝶优化QSAR特征选择

4.1 二进制帝王蝶特征选择算法

传统的帝王蝶优化算法的目标是在连续数据空间中寻找最优解.但特征选择是在二进制空间中搜索最优解.为了解决QSAR特征选择的问题，每个个体均表示为一个二进制字符串.字符串的长度等于特征的数量.字符串中的“1”表示选择了对应的特征，“0”表示未选择对应的特征.同时，偏最小二乘法用于计算每个个体的适应度值.

二进制帝王蝶优化(Binary Monarch Butterfly Optimization，BMBO)算法的过程与标准帝王蝶优化算法相似，只在两个地方做了修改：

1)在初始化阶段，所有个体的位置均随机生成，然后根据给定阈值转换成二进制字符串.

2)因在调整算子中获得的值是连续值，使用Tanh()传递函数[14]将其映射为二进制值.故当rand>BAR时，公式(5)修改为：

(7)

需要注意的是，为了适应最小值优化，本研究用准确度的负值即-Q2作为适应度函数.

4.2 快速非支配排序策略

为了同时满足最小化准确度负值和最小化特征数的要求，引入了非支配排序算法.

首先，快速非支配排序涉及到的概念有如下两点：

1)对于两个向量x和y，当且仅当：

(8)

我们称x支配y，表示为xy，否则，这两个向量互不支配.

2)对于一个解x∈X，当且仅当：

{∃x∈X|xjxi，i≠j}

(9)

则称该解x为帕累托最优.所有帕累托最优解的集合称为帕累托解集，定义如式(10)所示：

Ps={xi，xj∈X|∃xjxi，i≠j}

(10)

一个包含帕累托解集目标函数值的集合称为帕累托前沿，表示为式(11)：

Pf={f(xi)|xi∈Ps}

(11)

因此，快速非支配排序算法对每个解x均需计算两个参数：1)种群中支配x的个体的数目nx；2)种群中被x支配的个体的集合Sx.

快速非支配排序的整个过程描述如下：

1)找到nx=0的个体，然后将这些个体存储到当前集合Fi中，其中i代表级别；

2)找出在Fi中由x支配的个体的集合Sx.对于Sx中的每个个体x，用nx-1更新nx.如果nx=0，则将x存储在集合H中.

3)F1中的个体为第一非支配水平的个体.H代表当前集合.重复上述这些步骤，直到所有的个体被分配到不同级别的集合为止.

同时，还对调整算子进行了修改.与传统的非支配排序算法不同，我们只考虑第一层帕累托前沿F1，并且只计算F1中个体的拥挤度.当rand≤p时，调整算子中的个体按照式(12)进行更新：

(12)

4.3 准确度优先策略

该策略的思想是当子代个体拥有更好的精度时，对应的父代个体将会被其取代.增加该策略的原因是：在QSAR建模中，模型的准确度指标优先于特征数指标.采用准确度优先策略可以避免多目标算法耗费大量计算资源来搜索低准确度的空间，即图1中矩形框内的区域.

图1 低准确度的区域示意图Fig.1 Illustration of region with the low accuracy

4.4 突变策略

为了解决早熟的问题，整个种群分为子组1、子组2和子组3，每个子组的个体数量均为NP/3.子组1中的突变率设置为随着迭代次数的增加而降低，子组2的突变率则设置为随机数，子组3中的个体不进行突变.

4.5 准确度优先的多目标帝王蝶优化QSAR特征选择的基本框架

综合上述3个策略，准确度优先的多目标帝王蝶优化(Accuracy Prior Multi-Objective Monarch Butterfly Optimization，APMOMBO)QSAR特征选择的步骤如算法2所示.该算法最终输出的最优粒子是准确度最高的粒子.

算法2.准确度优先的多目标帝王蝶优化QSAR特征选择算法伪代码

初始化参数：种群大小NP，子种群1的个体占比p，调整率BAR，迁移周期peri，最大步长Smax，t=1，最大评价次数T；

随机初始化大小为NP的种群；

通过偏最小二乘法计算种群中每个个体的适应度值，并统计每个个体所选择的特征数量；

通过非支配排序策略对种群中所有个体进行排序；

Whilet≤Tdo

根据p将种群分成子种群1和子种群2；

For子种群1中的每个个体

根据公式(2)进行个体的更新；

End

For子种群2中的每个个体

根据公式(7)和公式(12)进行个体的更新；

End

将两个更新后的子种群合并成一个种群，并将该种群分成子组1、子组2和子组3；

根据突变策略对前两个子组中的个体进行变异操作；

通过偏最小二乘法计算种群中每个个体的适应度值，并统计每个个体所选择的特征数量；

根据准确度优先策略更新子组1和子组2中的父代个体；

将所有子组合并为一个种群；

通过非支配排序策略对种群中所有个体进行排序；

t=t+1；

End

输出最优个体的位置xbest及其对应的适应度值f(xbest).

5 算法性能测试

5.1 实验设置及性能指标

为了测试我们提出的算法的特征选择性能，本文选择了文献[11]中的3个基准数据集，如表1所示.

表1 测试数据集Table 1 List of datasets used in the experiments

同时，将提出的APMOMBO算法和BMBO、PSO和WS-PSO进行比较.在本文中BMBO和APMOMBO算法的相关参数设置为：NP=150，p=0.5，BAR=0.5，peri=1.2，Smax=1.0，T=20000.除种群数量和评价次数分别为150和20000外，PSO和WS-PSO算法的其它参数设置与文献[11]相同，即PSO的学习率a=0.5，WS-PSO的学习率a=0.5，b=0.8，突变率c=0.05以及加权抽样比例为α=0.5.偏最小二乘法中潜在化合物的最大数量被设置为20，并且通过5折交叉验证来确定潜在化合物的最佳数量.为了消除随机性，所有算法在3个数据集上均独立运行100次.本文所有实验在硬件配置为Intel(R)Core(TM)i5-7500 CPU @3.40GHz、内存为16.0GB的计算机上进行，操作系统为Microsoft Windows 10.编程语言为MATLAB，编译环境为MATLAB R2018b.我们使用文献[11]中描述的5折交叉验证方法来评估QSAR特征选择的性能.本文的性能指标包括准确度Q2、均方根误差[15](表示为RMSECV)和特征数量(表示为NSD).准确度的值越大越好，其它两个指标则越小越好.同时，我们还计算了各指标的标准差.

5.2 性能比较与分析

由于BMBO未在QSAR建模特征选择中得到应用，我们给出了BMBO和PSO和WS-PSO的比较结果，如表2所示，其中粗体为较优结果.在准确度Q2和均方根误差RMSECV方面，BMBO的性能均优于PSO及WS-PSO算法，且标准差相对更低.这表明将BMBO算法用于QSAR建模特征选择，可以保持较高的准确度和稳定性.但BMBO算法在Artemisinin数据集和Selwood数据集的降维效果较差，这说明需要进一步改进BMBO的降维性能.

为了评估不同策略与BMBO算法结合的性能和效率，我们将3个策略逐个添加到BMBO算法中，得到3个不同的算法，即仅采用快速非支配排序策略的二进制帝王蝶优化算法(称为IBMBO-1)，同时采用快速非支配排序策略和准确度优先策略的二进制帝王蝶优化算法(称为IBMBO2)，以及同时采用快速非支配排序策略、准确度优先策略以及突变策略的二进制帝王蝶优化算法即APMOMBO.表3给出了BMBO和上述不同策略结合的实验结果，其中粗体为较优结果.

从表3可以看出，IBMBO-2和APMOMBO在准确度和均方根误差方面比IBMBO-1表现更好，而在特征数量方面比IBMBO-1表现得差.这表明，准确度优先策略和准确度优先+突变策略有利于提高特征选择的准确性和稳定性，仅使用快速非支配排序算法很难搜索到所有有用的特征来保持较高的准确性.与IBMBO-2相比，APMOMBO除了在Selwood数据集的特征数量之外的所有性能指标上都有提升.结果表明，采用突变策略避免了早熟，提高了精度，减少了均方根误差和特征数.

对比表2和表3中的数据可以看出，APMOMBO在3个性能指标上都比BMBO、WS-PSO和PSO表现得更好.这说明改进的BMBO算法通过使用快速非支配排序、准确度优先和突变策略，提高了精度，减少了均方根误差和特征数量.

表2 3种算法用于QASR建模特征选择的性能指标对比Table 2 Comparison of the algorithms for feature selection in QSAR

表3 BMBO算法结合3种策略的实验结果Table 3 Experimental results of BMBO algorithms with three strategies

6 总 结

在本文中，我们应用二进制帝王蝶优化来选择分子描述符用于QSAR建模，并且结合快速非支配排序、准确度优先和突变3种策略来提高BMBO的性能.

对BMBO和准确度优先的多目标帝王蝶优化QSAR特征选择方法在3个基准数据集上进行了测试.实验结果表明BMBO算法可以用于QSAR特征选择，其在准确度和均方根误差上的性能均优于PSO和WS-PSO，但在降维方面略有不足.而准确度优先的多目标帝王蝶优化QSAR特征选择方法则在所有指标上均优于PSO、WS-PSO和BMBO算法.这表明快速非支配排序、准确度优先和突变策略提高了特征选择的准确度，降低了均方根误差，并减少了特征的数量.该算法可有效实现QSAR建模中的特征选择.