薛艷鋒, 劉繼華, 張 翔, 薛志文

(呂梁學(xué)院計算機(jī)科學(xué)與技術(shù)系, 山西 呂梁 033000)

1 引言(Introduction)

關(guān)于流行病動力學(xué)的數(shù)學(xué)模型是了解和干預(yù)流行病在現(xiàn)實世界傳播的必要工具[1]。比如,SIS(Susceptible-Infected-Susceptible)模型可以刻畫個體被感染、恢復(fù)再被感染、恢復(fù)的周而復(fù)始的流行病傳播過程[2]。接觸網(wǎng)絡(luò)(節(jié)點表示個體,連邊表示個體之間的接觸)可以更真實地刻畫現(xiàn)實世界流行病的傳播過程[3]。在SIS模型中,如果感染強(qiáng)度大于傳染閾值,則傳染病會流行,反之,則不會流行;其中,傳染閾值為接觸網(wǎng)絡(luò)譜半徑的倒數(shù)[4]。然而現(xiàn)實世界中,不可能存在一個孤立的群體(用接觸網(wǎng)絡(luò)G表示)。如果接觸網(wǎng)絡(luò)G1與另一個接觸網(wǎng)絡(luò)G2相互連接,則彼此的傳染閾值都會減小[4],進(jìn)而可能導(dǎo)致相互連接之前都不會暴發(fā)的流行病在連接之后暴發(fā)。所以,在連接之后仍然不會暴發(fā)流行病的條件下,求取最大連邊數(shù)很有必要。該過程可通過梯度下降算法計算,但是矩陣到譜半徑的映射無法用可微函數(shù)表示,所以無法直接應(yīng)用。本文利用目標(biāo)矩陣的F范數(shù)(該F范數(shù)大于譜半徑且是矩陣的凸函數(shù))構(gòu)建損失函數(shù),并初始化連接矩陣的參數(shù)滿足特定分布,進(jìn)而通過梯度下降算法成功找到滿足約束條件的最優(yōu)解。

2 預(yù)備知識與問題定義(Preliminary knowledge and definition of problem)

2.1 基于接觸網(wǎng)絡(luò)的SIS模型

在基于單個接觸網(wǎng)絡(luò)G(用鄰接矩陣A表示,元素為1表示接觸,為0表示未接觸)的SIS模型中,個體分為感染者(I)和易感者(S),如圖1所示,恢復(fù)率δ為感染者恢復(fù)為易感者的概率(如感染者2、4所示),感染率β為感染者通過接觸(用連邊表示)成功感染鄰居易感者的概率(如感染者7成功感染鄰居易感者4)。其中,感染強(qiáng)度τ和傳染閾值ε定義如下:

圖1 基于接觸網(wǎng)絡(luò)的SIS模型Fig.1 SIS models based on the contact networks

τ=β/δ

(1)

ε=1/ρ(A)

(2)

其中,ρ(A)為接觸網(wǎng)絡(luò)G的譜半徑。

如圖2所示,當(dāng)接觸網(wǎng)絡(luò)G1(黑色表示)與另一個接觸網(wǎng)絡(luò)G2(灰色表示)相互連接時,分別用A11和A22表示接觸網(wǎng)絡(luò)G1和G2的鄰接矩陣,N1和N2表示節(jié)點個數(shù),δ1和δ2表示恢復(fù)率,A12表示G1和G2之間的個體接觸(黑色虛線表示)。βij表示Gj對Gi的感染率(i,j∈{1,2})。此外,有以下幾點需要補(bǔ)充。

圖2 相互連接的接觸網(wǎng)絡(luò)Fig.2 An interconnected contact network

(2)由于感染的過程同時涉及感染者和易感者,所以感染率β可分解為感染者呼出病毒的速率μ和易感者吸入病毒的概率ω的乘積,則有

β11=μ1ω1

(3)

β22=μ2ω2

(4)

2.2 定義

無標(biāo)度網(wǎng)絡(luò)[5]:通過生長與偏好連接生成的網(wǎng)絡(luò),具體過程為每次輸入m條邊與已經(jīng)存在的m個節(jié)點相連(生長),偏好連接體現(xiàn)在優(yōu)先連接度大的節(jié)點。由于無標(biāo)度網(wǎng)絡(luò)更能代表真實網(wǎng)絡(luò)[6],所以本文設(shè)定G1和G2都為無標(biāo)度網(wǎng)絡(luò),通過networkx庫[7]生成。

譜半徑:矩陣A∈N×N的譜半徑表達(dá)式為ρ(A)=max(|λ1|,…,|λN|),即矩陣特征值模的最大值,又因為本文假定所有網(wǎng)絡(luò)為無向網(wǎng)絡(luò),即鄰接矩陣為對稱矩陣,所以譜半徑可簡化為鄰接矩陣的最大特征值。

矩陣的F范數(shù)‖A‖F(xiàn):把矩陣中每個元素的平方求和再開根號,即

(5)

2.3 問題

G1和G2相互連接之前,感染強(qiáng)度和傳染閾值分別如下:

τ11=β11/δ1

(6)

τ22=β22/δ2

(7)

ε11=1/ρ(A11)

(8)

ε22=1/ρ(A22)

(9)

其中,當(dāng)τ11<ε11且τ22<ε22時,G1和G2都不會暴發(fā)流行病。文獻(xiàn)[4]指出:G1連接G2之后,G1的傳染閾值下降:

ε11,c=1/ρ(H)

(10)

其中,

(11)

(12)

s.t.ε11,c>τ11

(13)

依據(jù)公式(10)可轉(zhuǎn)化為如下最優(yōu)化問題:

(14)

s.t.ρ(H)<1/τ11

(15)

3 算法(Algorithm)

深度學(xué)習(xí)的函數(shù)擬合過程都是利用損失函數(shù)通過梯度下降算法調(diào)整參數(shù)矩陣,從而在參數(shù)空間找到一組最優(yōu)參數(shù)的過程[8-11]。同理,本文也是利用損失函數(shù)在兩個接觸網(wǎng)絡(luò)G1和G2的接觸連邊A12空間中尋找一組最多連邊的組合。

設(shè)矩陣H的譜半徑和F范數(shù)與G1感染強(qiáng)度的倒數(shù)差值分別如下:

dρ(H)=ρ(H)-1/τ11

(16)

d‖H‖F(xiàn)=‖H‖F(xiàn)-1/τ11

(17)

由于矩陣到譜半徑的映射無法用可微函數(shù)計算,所以本文采用公式(17)利用F范數(shù)替代公式(16)構(gòu)建損失函數(shù),具體流程如下。

(1)初始化矩陣WN1×N2的元素滿足區(qū)間為[1,2]的均勻分布,替換矩陣A12代入公式(11)可得矩陣HW。

(2)設(shè)置矩陣A12=W>0.5,即

(18)

其中,i∈{1,2,…,N1},j∈{1,2,…,N2}。代入公式(11)可得矩陣H。此時,由于矩陣W>A12(矩陣A12表示元素全為1的矩陣),所以可得矩陣不等式H

(3)設(shè)計初始化參數(shù)矩陣WN1×N2的指導(dǎo)原則是使G1的感染強(qiáng)度倒數(shù)小于矩陣H的譜半徑ρ(H),即有

1/τ11<ρ(H)

(19)

進(jìn)而,不等式

0<1/τ11<ρ(H)<‖HW‖F(xiàn)

(20)

成立。

4 實驗結(jié)果及分析(Experimental results and analysis)

4.1 系數(shù)設(shè)置

本文選擇的系數(shù)設(shè)置如表1所示,N表示節(jié)點數(shù),m表示無標(biāo)度網(wǎng)絡(luò)每次生長的邊數(shù),τ22表示G2的感染強(qiáng)度,“—”表示占位符,不代表任何實際含義。

取兩種極端情況如下。

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感染強(qiáng)度的倒數(shù)1/τ11范圍:

1/τ11∈(1/maxε11,c,1/minε11,c)=(7.53,105.80)

(25)

基于此,本文選擇G1感染強(qiáng)度的倒數(shù)1/τ11分別為20、40、60、80。

表 1 系數(shù)設(shè)置Tab.1 Coefficients settings

4.2 實驗方案

4.3 實驗結(jié)果

如表2所示,在四種G1感染強(qiáng)度的倒數(shù)1/τ11的不同取值情況下,梯度下降的方法計算得到滿足條件的連邊數(shù)最大及標(biāo)準(zhǔn)差最小。同時,當(dāng)G1的感染強(qiáng)度越大(越容易暴發(fā)流行病,即1/τ11越小),本文所提算法取得的效果越明顯。

表 2 實驗結(jié)果對比Tab.2 Comparison of experimental results

4.4 理論分析

圖3 當(dāng)1/τ11=20時,譜半徑和F范數(shù)與感染強(qiáng)度倒數(shù)變化示意圖Fig.3 Schematic diagram of reciprocal changes of spectral radius and Fnorm with infection intensity when 1/τ11=20

4.5 算法運(yùn)行過程描述

從圖3可以看出,代表H矩陣譜半徑ρ(H)的曲線在訓(xùn)練初期,譜半徑的值處于最大值105.8,此時對應(yīng)的連邊矩陣A12的元素都為1,即G1與G2中的每個個體之間都相互連接。隨著訓(xùn)練的持續(xù)進(jìn)行,譜半徑ρ(H)開始下降,即去邊的過程開始(開始時的坐標(biāo)通過元組顯示),直到優(yōu)化過程結(jié)束。

4.6 算法的擴(kuò)展性分析

本算法的擴(kuò)展性從三個方面進(jìn)行分析。

(1)算法層面:如果A12表示接觸概率,則其他三種連邊方式失效,本文算法只需要把公式(18)替換為(A12)ij=sigmoid(Wij)即可。

(2)硬件層面:可以充分利用圖形處理器(Graphics Processing Unit,GPU)并行加速計算,對于規(guī)模較大的網(wǎng)絡(luò),可以明顯提高運(yùn)算效率。

(3)軟件層面:各種深度學(xué)習(xí)的開源框架TensorFlow[12]和PyTorch[13]不僅支持GPU運(yùn)算,而且提供了一系列完整的損失函數(shù)、初始化參數(shù)矩陣、自適應(yīng)學(xué)習(xí)率算法等,這些都為本文所提算法的擴(kuò)展提供了便利。

5 結(jié)論(Conclusion)

本文利用矩陣的F范數(shù)替代譜半徑構(gòu)建損失函數(shù),成功地解決了兩個接觸網(wǎng)絡(luò)的最多連邊優(yōu)化問題,并從實驗結(jié)果和理論分析證明了所提算法的有效性。而且,通過算法運(yùn)行過程的描述,直觀地展示了尋找最多連邊的整個過程并能對關(guān)鍵環(huán)節(jié)給出解釋。此外,從算法層面、硬件層面以及軟件層面對所提算法的擴(kuò)展性進(jìn)行了分析,為解決不可微損失函數(shù)的優(yōu)化問題提供了新的思路。