薛艳锋, 刘继华, 张 翔, 薛志文
(吕梁学院计算机科学与技术系, 山西 吕梁 033000)
关于流行病动力学的数学模型是了解和干预流行病在现实世界传播的必要工具[1]。比如,SIS(Susceptible-Infected-Susceptible)模型可以刻画个体被感染、恢复再被感染、恢复的周而复始的流行病传播过程[2]。接触网络(节点表示个体,连边表示个体之间的接触)可以更真实地刻画现实世界流行病的传播过程[3]。在SIS模型中,如果感染强度大于传染阈值,则传染病会流行,反之,则不会流行;其中,传染阈值为接触网络谱半径的倒数[4]。然而现实世界中,不可能存在一个孤立的群体(用接触网络G表示)。如果接触网络G1与另一个接触网络G2相互连接,则彼此的传染阈值都会减小[4],进而可能导致相互连接之前都不会暴发的流行病在连接之后暴发。所以,在连接之后仍然不会暴发流行病的条件下,求取最大连边数很有必要。该过程可通过梯度下降算法计算,但是矩阵到谱半径的映射无法用可微函数表示,所以无法直接应用。本文利用目标矩阵的F范数(该F范数大于谱半径且是矩阵的凸函数)构建损失函数,并初始化连接矩阵的参数满足特定分布,进而通过梯度下降算法成功找到满足约束条件的最优解。
2.1 基于接触网络的SIS模型
在基于单个接触网络G(用邻接矩阵A表示,元素为1表示接触,为0表示未接触)的SIS模型中,个体分为感染者(I)和易感者(S),如图1所示,恢复率δ为感染者恢复为易感者的概率(如感染者2、4所示),感染率β为感染者通过接触(用连边表示)成功感染邻居易感者的概率(如感染者7成功感染邻居易感者4)。其中,感染强度τ和传染阈值ε定义如下:
图1 基于接触网络的SIS模型Fig.1 SIS models based on the contact networks
τ=β/δ
(1)
ε=1/ρ(A)
(2)
其中,ρ(A)为接触网络G的谱半径。
如图2所示,当接触网络G1(黑色表示)与另一个接触网络G2(灰色表示)相互连接时,分别用A11和A22表示接触网络G1和G2的邻接矩阵,N1和N2表示节点个数,δ1和δ2表示恢复率,A12表示G1和G2之间的个体接触(黑色虚线表示)。βij表示Gj对Gi的感染率(i,j∈{1,2})。此外,有以下几点需要补充。
图2 相互连接的接触网络Fig.2 An interconnected contact network
(2)由于感染的过程同时涉及感染者和易感者,所以感染率β可分解为感染者呼出病毒的速率μ和易感者吸入病毒的概率ω的乘积,则有
β11=μ1ω1
(3)
β22=μ2ω2
(4)
2.2 定义
无标度网络[5]:通过生长与偏好连接生成的网络,具体过程为每次输入m条边与已经存在的m个节点相连(生长),偏好连接体现在优先连接度大的节点。由于无标度网络更能代表真实网络[6],所以本文设定G1和G2都为无标度网络,通过networkx库[7]生成。
谱半径:矩阵A∈N×N的谱半径表达式为ρ(A)=max(|λ1|,…,|λN|),即矩阵特征值模的最大值,又因为本文假定所有网络为无向网络,即邻接矩阵为对称矩阵,所以谱半径可简化为邻接矩阵的最大特征值。
矩阵的F范数‖A‖F:把矩阵中每个元素的平方求和再开根号,即
(5)
2.3 问题
G1和G2相互连接之前,感染强度和传染阈值分别如下:
τ11=β11/δ1
(6)
τ22=β22/δ2
(7)
ε11=1/ρ(A11)
(8)
ε22=1/ρ(A22)
(9)
其中,当τ11<ε11且τ22<ε22时,G1和G2都不会暴发流行病。文献[4]指出:G1连接G2之后,G1的传染阈值下降:
ε11,c=1/ρ(H)
(10)
其中,
(11)
(12)
s.t.ε11,c>τ11
(13)
依据公式(10)可转化为如下最优化问题:
(14)
s.t.ρ(H)<1/τ11
(15)
深度学习的函数拟合过程都是利用损失函数通过梯度下降算法调整参数矩阵,从而在参数空间找到一组最优参数的过程[8-11]。同理,本文也是利用损失函数在两个接触网络G1和G2的接触连边A12空间中寻找一组最多连边的组合。
设矩阵H的谱半径和F范数与G1感染强度的倒数差值分别如下:
dρ(H)=ρ(H)-1/τ11
(16)
d‖H‖F=‖H‖F-1/τ11
(17)
由于矩阵到谱半径的映射无法用可微函数计算,所以本文采用公式(17)利用F范数替代公式(16)构建损失函数,具体流程如下。
(1)初始化矩阵WN1×N2的元素满足区间为[1,2]的均匀分布,替换矩阵A12代入公式(11)可得矩阵HW。
(2)设置矩阵A12=W>0.5,即
(18)
其中,i∈{1,2,…,N1},j∈{1,2,…,N2}。代入公式(11)可得矩阵H。此时,由于矩阵W>A12(矩阵A12表示元素全为1的矩阵),所以可得矩阵不等式H (3)设计初始化参数矩阵WN1×N2的指导原则是使G1的感染强度倒数小于矩阵H的谱半径ρ(H),即有 1/τ11<ρ(H) (19) 进而,不等式 0<1/τ11<ρ(H)<‖HW‖F (20) 成立。 本文选择的系数设置如表1所示,N表示节点数,m表示无标度网络每次生长的边数,τ22表示G2的感染强度,“—”表示占位符,不代表任何实际含义。 取两种极端情况如下。 minρ(H)=ρ(G1)=7.53 (21) 这种情况下,G1的传染阈值最大,即 maxε11,c=1/7.53 (22) 此时,表示G1与G2没有任何连接,即针对矩阵H的研究等价于研究完全孤立的群体G1,所以这种情况不属于本文研究的范围。 maxρ(H)=ρ(G1)=105.80 (23) 这种情况下,G1的传染阈值下降为最小值: minε11,c=1/105.80 (24) 如果在这种情况下仍然满足前提条件τ11 综合上述两种情况,本文研究的G1感染强度的倒数1/τ11范围: 1/τ11∈(1/maxε11,c,1/minε11,c)=(7.53,105.80) (25) 基于此,本文选择G1感染强度的倒数1/τ11分别为20、40、60、80。 表 1 系数设置Tab.1 Coefficients settings 如表2所示,在四种G1感染强度的倒数1/τ11的不同取值情况下,梯度下降的方法计算得到满足条件的连边数最大及标准差最小。同时,当G1的感染强度越大(越容易暴发流行病,即1/τ11越小),本文所提算法取得的效果越明显。 表 2 实验结果对比Tab.2 Comparison of experimental results 图3 当1/τ11=20时,谱半径和F范数与感染强度倒数变化示意图Fig.3 Schematic diagram of reciprocal changes of spectral radius and Fnorm with infection intensity when 1/τ11=20 从图3可以看出,代表H矩阵谱半径ρ(H)的曲线在训练初期,谱半径的值处于最大值105.8,此时对应的连边矩阵A12的元素都为1,即G1与G2中的每个个体之间都相互连接。随着训练的持续进行,谱半径ρ(H)开始下降,即去边的过程开始(开始时的坐标通过元组显示),直到优化过程结束。 本算法的扩展性从三个方面进行分析。 (1)算法层面:如果A12表示接触概率,则其他三种连边方式失效,本文算法只需要把公式(18)替换为(A12)ij=sigmoid(Wij)即可。 (2)硬件层面:可以充分利用图形处理器(Graphics Processing Unit,GPU)并行加速计算,对于规模较大的网络,可以明显提高运算效率。 (3)软件层面:各种深度学习的开源框架TensorFlow[12]和PyTorch[13]不仅支持GPU运算,而且提供了一系列完整的损失函数、初始化参数矩阵、自适应学习率算法等,这些都为本文所提算法的扩展提供了便利。 本文利用矩阵的F范数替代谱半径构建损失函数,成功地解决了两个接触网络的最多连边优化问题,并从实验结果和理论分析证明了所提算法的有效性。而且,通过算法运行过程的描述,直观地展示了寻找最多连边的整个过程并能对关键环节给出解释。此外,从算法层面、硬件层面以及软件层面对所提算法的扩展性进行了分析,为解决不可微损失函数的优化问题提供了新的思路。
4.1 系数设置
4.2 实验方案
4.3 实验结果
4.4 理论分析
4.5 算法运行过程描述
4.6 算法的扩展性分析