朱 霞，潘瑩麗

（1.中南財(cái)經(jīng)政法大學(xué) 統(tǒng)計(jì)與數(shù)學(xué)學(xué)院，武漢 430073；2.湖北大學(xué)a.數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院；b.應(yīng)用數(shù)學(xué)湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430062）

0 引言

分布式計(jì)算的主要特點(diǎn)是對海量數(shù)據(jù)進(jìn)行分布式數(shù)據(jù)挖掘。在分布式系統(tǒng)中，不同的機(jī)器之間需要數(shù)據(jù)的傳輸。在主從模式的分布式系統(tǒng)中，研究者普遍關(guān)注且不可避免的一個(gè)問題是由Master 和Worker 之間信息傳輸所導(dǎo)致的溝通成本問題。為了降低溝通成本，一種常用的方法是Zhang 等（2013）[1]提出的one-shot 方法。該方法讓每個(gè)Worker先根據(jù)自己的數(shù)據(jù)得到一個(gè)比較好的估計(jì)值，再把每個(gè)估計(jì)值由不同的Worker傳遞給Master，然后Master將各個(gè)Worker 傳遞過來的估計(jì)值求平均得最到終的估計(jì)量。one-shot 方法只需進(jìn)行一輪主從機(jī)器之間信息的傳輸，且傳輸?shù)膬H是一個(gè)向量。然而，該方法為了使平均估計(jì)量達(dá)到一定的收斂速度，要求每個(gè)機(jī)器必須獲得至少個(gè)樣本，機(jī)器的數(shù)量必須比小得多[2,3]。Jordan 等（2019）[4]提出的基于交互有效的替代損失（CSL）函數(shù)的分布式學(xué)習(xí)方法較好地克服了one-shot 方法的缺點(diǎn)。在頻域框架中，CSL函數(shù)可以像全局損失函數(shù)一樣用于形成點(diǎn)估計(jì)或進(jìn)行統(tǒng)計(jì)推斷，且在有限次的溝通次數(shù)下，基于CSL 函數(shù)的估計(jì)量與基于全局損失函數(shù)的估計(jì)量具有相同的收斂速度。借鑒構(gòu)造CSL函數(shù)的思想，Pan（2021）[5]和潘瑩麗等（2022）[6]分別基于大規(guī)模高維尾期望回歸模型和大規(guī)模Huber 回歸模型給出了一個(gè)溝通有效的分布式優(yōu)化方法。然而，這兩篇文章給出的解決方案只是每個(gè)Worker 僅計(jì)算梯度，由Master 求解優(yōu)化問題，從而導(dǎo)致Master 一直在工作，而其他Worker 在空轉(zhuǎn)。因此，在主從分布式架構(gòu)下，如何避免將任務(wù)集中于Master 一臺(tái)機(jī)器，進(jìn)而研究如何改進(jìn)分布式優(yōu)化方法非常必要。

受交互有效的替代損失（CSL）函數(shù)的啟發(fā)[4]，本文基于尾期望回歸模型提出一種廣義替代損失函數(shù)——梯度增強(qiáng)型損失（GEL）函數(shù)，能有效避免CSL函數(shù)導(dǎo)致的Worker空轉(zhuǎn)而Master一直工作的局面。為了對參數(shù)進(jìn)行估計(jì)和對變量進(jìn)行選擇，本文對GEL 函數(shù)施加平滑裁剪絕對偏差（SCAD）懲罰和自適應(yīng)最小絕對收縮與選擇算子（adaptive LASSO）懲罰[7,8]，并構(gòu)造正則化梯度增強(qiáng)型損失（GEL）函數(shù)。由于Boyd等（2011）[9]指出交替方向乘子法（ADMM）算法適用于大規(guī)模統(tǒng)計(jì)推斷和分布式凸優(yōu)化問題，故為了對正則化GEL函數(shù)進(jìn)行優(yōu)化，本文還提出一種增廣型極大交替方向乘子（Proximal-ADMM）算法。本文提出的方法運(yùn)用正則化GEL函數(shù)，在計(jì)算上充分利用每臺(tái)機(jī)器（Worker）的計(jì)算優(yōu)勢，加快收斂速度，使得所有機(jī)器（Worker）可以并行優(yōu)化各自對應(yīng)的GEL函數(shù)，然后傳輸?shù)街鳈C(jī)（Master）上對結(jié)果進(jìn)行整合得到最終的估計(jì)。在溝通成本方面，本文將通過模擬和實(shí)證研究驗(yàn)證Proximal-ADMM 算法可以在有限輪Master 和Worker 之間的溝通中收斂于基于全樣本（所有數(shù)據(jù)）的Centralize方法的結(jié)果。

1 高維尾期望回歸模型分布式估計(jì)的改進(jìn)過程

1.1 高維尾期望回歸模型

其中，xi是p維協(xié)變量；β0(τ)是p維參數(shù)；εi(τ)是隨機(jī)誤差項(xiàng)，滿足給定xi，εi(τ)的第τ個(gè)尾期望為零。為了便于表示，以下將β0(τ)和εi(τ)分別簡記為β0和εi。假設(shè)參數(shù)真值β0是稀疏的，即若記Γ={k:βk0≠0} ，其中βk0是β0的第k個(gè)元素，q=|Γ|是集合Γ 中元素的個(gè)數(shù)，則稀疏假設(shè)意味著q?p。在大數(shù)據(jù)背景下，受單機(jī)存儲(chǔ)性能的制約，很難將全樣本存儲(chǔ)于一臺(tái)機(jī)器之中。在主從分布式架構(gòu)下，將數(shù)據(jù)隨機(jī)均勻存儲(chǔ)是解決單臺(tái)機(jī)器存儲(chǔ)瓶頸的重要方法。不妨設(shè)樣本總量N=nm，即有m臺(tái)機(jī)器，每臺(tái)機(jī)器隨機(jī)存儲(chǔ)n個(gè)樣本。記第j臺(tái)機(jī)器上的觀測數(shù)據(jù)為，則所有樣本可被記為：

定義基于全部觀測數(shù)據(jù)的全局損失函數(shù)如下：

其中，ρτ(u)=τu2I(u＞0 )+(1 -τ)u2I(u≤0 )是尾期望回歸中的凸損失函數(shù)。類似地，定義基于存儲(chǔ)在第j臺(tái)機(jī)器上數(shù)據(jù)的局部損失函數(shù)如下：

在大數(shù)據(jù)背景下，為了同時(shí)達(dá)到參數(shù)估計(jì)和變量選擇的目的，若直接最小化正則化的全局損失函數(shù)QN(β)將產(chǎn)生較高的計(jì)算成本，且對單機(jī)的性能要求較高。因此，尋找全局損失函數(shù)的一個(gè)替代損失函數(shù)，將是解決單機(jī)存儲(chǔ)瓶頸和單機(jī)計(jì)算耗時(shí)耗力的有效方法之一。

1.2 正則化梯度增強(qiáng)型損失（GEL）函數(shù)

為了構(gòu)造全局損失函數(shù)的一個(gè)替代損失函數(shù)，將QN(β)在任意初始估計(jì)量處進(jìn)行一階泰勒展開得：

以下將替代損失函數(shù)Q(β)稱為GEL函數(shù)。

Jordan等（2019）[4]將R(β)替換成基于第1臺(tái)機(jī)器上數(shù)據(jù)的R1(β)，并提出交互有效的替代損失（CSL）函數(shù)如下：

由式（6）和式（7）可知，GEL 函數(shù)更廣泛，CSL 函數(shù)是GEL 函數(shù)的一個(gè)特例。在高維情形下，為避免過度擬合，提高模型的泛化能力，通常對GEL 函數(shù)施加一個(gè)懲罰項(xiàng)得到正則化GEL函數(shù)如下：

其中，pλ(·)是一個(gè)懲罰函數(shù)，λ是一個(gè)調(diào)優(yōu)參數(shù)。以下考慮兩種常用的懲罰：SCAD 懲罰和adaptive LASSO 懲罰。

對于SCAD懲罰，其懲罰函數(shù)形式為：

其中，a＞2，λ＞0。Fan和Li（2001）[7]建議取a=3.7，可通過交叉驗(yàn)證的方法選擇最優(yōu)參數(shù)λ。Zou 和Li（2008）[10]也提出SCAD懲罰函數(shù)可以做如下局部線性近似：

其 中，?pλ( ×) 是pλ(·) 的 一階導(dǎo)數(shù)，定義?pλ(u)=，符號e+表示若e＞0 ，則e+為e，否則為零。結(jié)合式（8）和SCAD 懲罰函數(shù)的局部線性近似式（9），省略與參數(shù)無關(guān)的常數(shù)項(xiàng)，構(gòu)建施加SCAD懲罰的GEL函數(shù)如下:

基于式（10），給出β的第一個(gè)估計(jì)量如下：

此外，本文將考慮adaptive LASSO 懲罰下參數(shù)的估計(jì)和變量選擇問題。將式（8）中的懲罰項(xiàng)換成adaptive LASSO懲罰[8]，得正則化GEL函數(shù)如下：

1.3 增廣型極大交替方向乘子（Proximal-ADMM）算法

為對正則化替代損失函數(shù)式（14）進(jìn)行優(yōu)化，現(xiàn)構(gòu)建增廣型Proximal-ADMM算法如下。

ADMM 算法通過將式（15）表示為以下的等效形式來解決：

其中，γ＞0 是調(diào)優(yōu)參數(shù)，為增廣項(xiàng)。根據(jù)凸優(yōu)化理論，優(yōu)化問題式（16）的拉格朗日函數(shù)為：

基于第j(j=1，2，…，m)臺(tái)機(jī)器（Worker）上的觀察值的優(yōu)化目標(biāo)函數(shù)L(β，z，θ)，ADMM 算法的迭代更新公式為：

其中，(βt+1，j，zt+1，j，θt+1，j)表示基于第j臺(tái)機(jī)器（Worker）上的數(shù)據(jù)計(jì)算出的第t+1 次迭代值。省略與參數(shù)無關(guān)的常數(shù)項(xiàng)，式（18）可表述為：

式（19）中的z步存在一個(gè)顯示解，且可以對z中的每個(gè)元素進(jìn)行優(yōu)化。事實(shí)上，對于i=1，2，…n，有：

其中，Proxρτ/nγ(·)為凸函數(shù)的近似點(diǎn)算子，定義如下：

對m臺(tái)機(jī)器（Worker）上的更新值zit+1，j求平均得z步的第t+1次更新值為：

對于一般的設(shè)計(jì)矩陣，β步的更新沒有顯示解。對式（19）中的β步增加一項(xiàng)，得到如下修正的β步：

其中，V=γ(ηIp×p-xTx)，η≥Λmax(xTx)，Ip×p是p×p單位矩陣，Λmax(·) 表示對實(shí)對稱矩陣取最大特征值，是由V定義的半內(nèi)積誘導(dǎo)的半范數(shù)。通過計(jì)算，可得：

由式（22），修正后的β步可變?yōu)椋?/p>

其中，j=1，2，…，m，x(k)表示設(shè)計(jì)矩陣x的第k(k=1，2，…，p)列。定義軟算子Shrink[u，a]=sgn(u)(|u|-a)+，其中sgn(·) 是符號函數(shù)，通過簡單計(jì)算，式（23）變?yōu)椋?/p>

對m臺(tái)機(jī)器（Worker）上的更新值βt+1，j求平均得β步的第t+1次更新值為：

同樣地，對基于m臺(tái)機(jī)器的式（19）中的θt+1，j取平均得θ步的第t+1次更新值為：

2 數(shù)值模擬

本文運(yùn)用Proximal-ADMM 算法（簡稱Paverage 算法）與基于第一臺(tái)機(jī)器所設(shè)計(jì)的分布式優(yōu)化方法（簡稱Psingle算法[5]）和基于全樣本執(zhí)行原始的ADMM 算法（簡稱Centralize算法）進(jìn)行模擬研究，來驗(yàn)證所提出的在SCAD懲罰和adaptive LASSO懲罰下尾期望回歸模型的增廣型Proximal-ADMM 算法的參數(shù)估計(jì)和變量選擇結(jié)果的有限樣本性能，以及其估計(jì)誤差與主從機(jī)器溝通次數(shù)的關(guān)系。

考慮回歸模型如下：

其中，協(xié)變量x1，x2，…，xp通過兩個(gè)步驟生成。步驟1：生成，協(xié)方差矩陣Σ的第(ij)元σij=0.5|i-j|，i，j=1，2，…，p；參數(shù)維數(shù)p=300。步驟2：設(shè)定，其中Φ(·)是標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。參數(shù)真值β6=β12=β15=β20=1，參數(shù)真值β的其余元素均取值為0。隨機(jī)誤差ε～N( 0，1)或ε～t(3)。樣本量N=2000，將其隨機(jī)均勻分配到m=20臺(tái)機(jī)器（Worker）上，每臺(tái)機(jī)器樣本量n=100?？紤]三個(gè)不同的尾期望值τ=0.3，0.5，0.7。將數(shù)據(jù)隨機(jī)分成訓(xùn)練集（Dtrain）和測試集（Dtest），采用20 折交叉驗(yàn)證法，通過優(yōu)化泛函來選 擇參 數(shù)λ，其中為基于訓(xùn)練集的參數(shù)估計(jì)量。

通過100次模擬，可從以下幾個(gè)方面比較上述Paverage、Centralize 和Psingle 三種方法的性能。分別給出SCAD懲罰和adaptive LASSO懲罰下的尾期望回歸模型的參數(shù)估計(jì)和變量選擇結(jié)果，如表1和表2所示。其中，Size為非零回歸系數(shù)的平均值≠0(k=1，2，…，p)。P1為在模擬時(shí)包含所有真正非零回歸系數(shù)的比例，即對于任何k≥1，估計(jì)值≠0且真實(shí)參數(shù)βk0≠0。P2為在模擬時(shí)選擇變量x1的比例。ER為估計(jì)誤差。在Size、P1、P2和ER的列中，括號內(nèi)的數(shù)字是該列指標(biāo)相應(yīng)的樣本標(biāo)準(zhǔn)差。

表1 SCAD懲罰下尾期望回歸的模擬結(jié)果

表2 adaptive LASSO懲罰下尾期望回歸的模擬結(jié)果

表1中的結(jié)果顯示，在參數(shù)估計(jì)方面，由Paverage方法和Psingle方法所生成的估計(jì)值可以與Centralize方法所生成的估計(jì)值相媲美。而Paverage方法的估計(jì)誤差比Psingle 方法的估計(jì)誤差要小。在變量選擇方面，當(dāng)τ=0.3 和τ=0.7 時(shí)，所有方法都成功地選取了5 個(gè)變量（x1、x6、x12、x15和x20），且有較高的選擇概率P1。對于協(xié)變量x1，除了τ=0.5 外，所有方法都展示出很高的選擇概率P2。表2中的結(jié)果同樣顯示了類似的結(jié)論。

為了證明所提出的Paverage、Centralize和Psingle三種方法是溝通有效的，本文模擬繪制出隨溝通次數(shù)變化的估計(jì)誤差趨勢圖，展示出施加SCAD懲罰的尾期望回歸隨溝通次數(shù)變化的估計(jì)誤差趨勢結(jié)果，如圖1 所示。施加adaptive LASSO懲罰的圖形和圖1類似。

圖1 SCAD懲罰下估計(jì)誤差與主從機(jī)器溝通次數(shù)的關(guān)系

圖1的結(jié)果表明，在經(jīng)過數(shù)次主從機(jī)器（Worker-Master）之間的溝通后，所提的Paverage 方法的估計(jì)誤差逐漸下降且趨近于Centralize 方法的估計(jì)誤差。此外，相對于Psingle方法的估計(jì)誤差，Paverage方法的估計(jì)誤差收斂到Centralize 方法的估計(jì)誤差的速度較快，表明基于GEL 函數(shù)的分布式優(yōu)化方法相對基于CSL 函數(shù)的分布式優(yōu)化方法在溝通有效方面具有更好的性能。

3 實(shí)證分析

本文將Paverage、Centralize和Psingle方法運(yùn)用于人類免疫缺陷病毒（HIV）藥物敏感性數(shù)據(jù)集的研究中，以驗(yàn)證所提出的方法在實(shí)際應(yīng)用中的可操作性。該數(shù)據(jù)集從網(wǎng)站http://hivdb.stanford.edu中獲取。研究人員指出，在未經(jīng)治療的艾滋病毒感染者中，每一個(gè)可能的單點(diǎn)突變每天發(fā)生104到105次，耐藥性已成為艾滋病毒治療的主要障礙。因此，了解突變對耐藥的影響是一個(gè)重要的研究課題。本文將主要對藥物efavirenz（EFV）敏感性數(shù)據(jù)進(jìn)行分析。在排除一些罕見突變后，該數(shù)據(jù)集包括1472 個(gè)HIV 分離株和197 個(gè)突變位點(diǎn)。將對藥物EFV 敏感性度量的自然對數(shù)作為響應(yīng)變量，將197個(gè)突變位點(diǎn)作為協(xié)變量。由于病毒突變的位置通常非常稀少，因此可通過施加懲罰項(xiàng)來獲得稀疏解。本文將建立施加懲罰的尾期望回歸模型研究197個(gè)突變位點(diǎn)對藥物EFV敏感性的影響。

采 用Paverage、Centralize 和Psingle 方 法 分 析HIV 數(shù)據(jù)并比較三種方法的性能?？紤]三個(gè)不同的尾期望值τ=0.3，0.5，0.7 。將整個(gè)數(shù)據(jù)集隨機(jī)分成兩部分：訓(xùn)練數(shù)據(jù)集（Dtrain）和測試數(shù)據(jù)集（Dtest），樣本量分別為1200和272。訓(xùn)練數(shù)據(jù)集在12 臺(tái)機(jī)器上隨機(jī)分割，每臺(tái)機(jī)器隨機(jī)存儲(chǔ)100 個(gè)樣本。基于訓(xùn)練集的1200 個(gè)樣本，在最小化損失函數(shù)的基礎(chǔ)上，采用12 折交叉驗(yàn)證法選擇調(diào)優(yōu)參數(shù)λ，然后利用測試集上的數(shù)據(jù)，通過PE=(1 272)計(jì)算預(yù)測誤差，其中(λ)是由訓(xùn)練集得到的參數(shù)估計(jì)。重復(fù)上述過程20次，表3匯總了所選重要解釋變量的預(yù)測誤差，此外，列括號里的數(shù)表示其相應(yīng)的樣本標(biāo)準(zhǔn)差。從預(yù)測誤差的角度來看，運(yùn)用Paverage方法得到的預(yù)測誤差均在合理范圍內(nèi)。

表3 用三種方法比較HIV藥物敏感性數(shù)據(jù)集性能的預(yù)測結(jié)果

將在τ=0.5 下通過20次隨機(jī)分割后各種方法選擇的重要變量的頻數(shù)匯總到表4 中。τ=0.3，0.7 的結(jié)果與τ=0.5 類似。從表4 可看出，除了個(gè)別基因，使用adaptive LASSO 懲罰的Paverage 方法選擇的變量結(jié)果與Centralize方法選擇的變量結(jié)果高度一致。總而言之，這些結(jié)果表明，在艾滋病毒藥物敏感性數(shù)據(jù)集研究中，本文提出的方法在變量選擇和預(yù)測性能兩個(gè)方面均有很好的效果。

表4 用三種方法比較HIV藥物敏感性數(shù)據(jù)集性能的頻率表（τ=0.5）

4 結(jié)論

在主從分布式架構(gòu)下，本文通過施加SCAD 懲罰和adaptive LASSO懲罰提出了施加懲罰的大規(guī)模尾期望回歸模型。為了降低機(jī)器之間的溝通成本，本文又構(gòu)造出一種新的全局損失函數(shù)——溝通有效的GEL函數(shù)。然后本文將基于全局損失函數(shù)的高維正則化估計(jì)問題轉(zhuǎn)化成了基于溝通有效的GEL函數(shù)的高維正則化估計(jì)問題。為了優(yōu)化施加懲罰的GEL 函數(shù)，本文還提出了增廣型Proximal-ADMM算法，在計(jì)算上，在Worker機(jī)器上計(jì)算估計(jì)量和在Master機(jī)器上整合估計(jì)量交替進(jìn)行，進(jìn)而加快了算法的收斂。在主從機(jī)器之間的溝通成本方面，經(jīng)過了幾輪機(jī)器之間的溝通后，本文提出的Paverage方法的估計(jì)誤差收斂于Centralize 方法的估計(jì)誤差，相對于Psingle 方法的估計(jì)誤差，Paverage方法的估計(jì)誤差收斂到Centralize方法的估計(jì)誤差的速度較快。數(shù)值模擬和實(shí)證分析的結(jié)果表明，本文提出的方法在樣本量有限的情況下具有良好的性能且在實(shí)際應(yīng)用中具有可操作性。