趙 磊 陳 玎 朱道立,

(1.上海交通大學(xué) 安泰經(jīng)濟(jì)與管理學(xué)院，上海 200030；2.上海交通大學(xué) 中美物流研究院，上海 200030)

(1)

即求出各特征(因素)的回歸系數(shù)構(gòu)成的向量β=(β1,…,βp)T∈Rp，和截距項(xiàng)β0。最經(jīng)典的回歸系數(shù)向量(β0,β)的估計(jì)模型為最小二乘估計(jì)模型(least squares estimator，LSE)：

(2)

隨著大數(shù)據(jù)技術(shù)的發(fā)展，數(shù)據(jù)產(chǎn)生的源頭逐漸增加，數(shù)據(jù)收集的難度逐漸降低，從而造成數(shù)據(jù)規(guī)模和復(fù)雜性逐漸增加。此外，數(shù)據(jù)的維度(屬性)通常很大，可以達(dá)到成百上千維甚至更高(例如文檔數(shù)據(jù)、基因表達(dá)數(shù)據(jù)、圖片數(shù)據(jù)、多媒體數(shù)據(jù)、貿(mào)易交易數(shù)據(jù)等)。在對應(yīng)的回歸問題中，樣本數(shù)據(jù)的維數(shù)p往往很大，樣本量n往往遠(yuǎn)遠(yuǎn)小于p，即p?n。以基因挖掘?yàn)槔?，通常是從成千上萬維的核苷酸數(shù)據(jù)中挖掘出臨床表現(xiàn)，但樣本量往往很小。例如Scheetz et.al.(2006)的例子中，需要在120只12周大的小白鼠的18976條基因信息中進(jìn)行挖掘，這一案例中p=18976，n=120。

在這種情況下，特征向量x_i的很多屬性往往線性相關(guān)，傳統(tǒng)的LSE方法通常很難處理這樣的數(shù)據(jù)。因此，一些統(tǒng)計(jì)學(xué)家在經(jīng)典LSE方法的基礎(chǔ)上加入變量篩選功能，以適應(yīng)高維數(shù)據(jù)的處理，即將LSE模型推廣為如下帶稀疏懲罰項(xiàng)的最小二乘形式：

(3)

其中，稀疏懲罰函數(shù)φλ(θ)可以有多種形式，常見的有嶺回歸(ridge regression)、套索回歸(least absolute shrinkage and selection operator，LASSO)、平滑削邊絕對偏離(smoothly clipped absolute deviation, SCAD)正則回歸等。隨著大數(shù)據(jù)時(shí)代的到來，這一類統(tǒng)計(jì)方法由于其在處理高維數(shù)據(jù)的優(yōu)越，越來越受到數(shù)據(jù)科學(xué)家的重視，近幾年在機(jī)器學(xué)習(xí)、統(tǒng)計(jì)學(xué)習(xí)、數(shù)據(jù)挖掘等領(lǐng)域被進(jìn)一步研究。

嶺回歸中，φλ(θ)通常選取如下平方形式：

φλ(θ)=λ|θ|2

(4)

LASSO回歸中，φλ(θ)通常選取如下絕對值形式：

φλ(θ)=λ|θ|

(5)

SCAD回歸中，φλ(θ)通常選取如下分段函數(shù)形式：

(6)

其中，a>2，λ>0。如圖1所示，二維空間中嶺回歸、LASSO和SCAD的懲罰函數(shù)分別為圖中虛線、點(diǎn)劃線、實(shí)線。

圖1 二維空間中嶺回歸、LASSO和SCAD的懲罰函數(shù)圖

前述三種帶稀疏懲罰項(xiàng)的回歸模型中，嶺回歸和LASSO回歸模型在統(tǒng)計(jì)性質(zhì)上是偏估計(jì)，SCAD回歸是近似無偏的，因此SCAD回歸具有更好的統(tǒng)計(jì)性質(zhì)，能夠被用來處理更加復(fù)雜的數(shù)據(jù)分析問題，例如存在奇異值的數(shù)據(jù)。但從回歸模型的求解上來講，LASSO回歸、嶺回歸相應(yīng)的最優(yōu)化模型是凸的，可用一階、二階算法，隨機(jī)坐標(biāo)下降算法等進(jìn)行求解(見Boyd)，而SCAD回歸是非凸的，造成了求解上的困難。但SCAD懲罰方法作為非凸懲罰項(xiàng)的前瞻性工作，直到現(xiàn)在仍受到理論和應(yīng)用領(lǐng)域的廣泛重視，因此有效求解SCAD的方法近年仍在被討論。在大數(shù)據(jù)背景下，待求解的回歸問題的數(shù)據(jù)量往往很大，數(shù)據(jù)源分布在不同的地理位置。這種情況使得進(jìn)行大數(shù)據(jù)分析時(shí)，往往難以獲得/處理完整數(shù)據(jù)。此外，由于數(shù)據(jù)規(guī)模大，算法的每一次迭代中更新全部變量會造成計(jì)算機(jī)的內(nèi)存不足。因此，在大數(shù)據(jù)環(huán)境下求解SCAD回歸問題時(shí)，往往難以實(shí)現(xiàn)在算法的每一次迭代更新全部變量。而現(xiàn)有求解SCAD回歸的算法(LQA、LLA、ADMM等)在每一次迭代時(shí)均需要更新全部變量，這造成了現(xiàn)有求解SCAD回歸的算法并不適用大數(shù)據(jù)背景下的應(yīng)用。隨機(jī)坐標(biāo)下降方法(stochastic coordinate descent, SCD)以其每一次迭代中僅需要更新一部分變量的流程設(shè)計(jì)，減少了計(jì)算的內(nèi)存需求，在大規(guī)模分布式最優(yōu)化問題中得到了廣泛的應(yīng)用。但目前理論上SCD算法只能處理帶凸懲罰項(xiàng)的回歸問題，由于SCAD回歸問題中的懲罰項(xiàng)是非凸和非光滑的，現(xiàn)有的隨機(jī)坐標(biāo)下降方法難以處理這一問題。

針對大數(shù)據(jù)時(shí)代高維數(shù)據(jù)、大規(guī)模數(shù)據(jù)、分布式存儲的特點(diǎn)，本文研究大規(guī)模SCAD回歸問題的隨機(jī)坐標(biāo)下降方法。為滿足算法收斂的要求，我們分析得到SCAD回歸模型的損失函數(shù)是導(dǎo)數(shù)Lipschitz、懲罰函數(shù)是semi-convex的分析結(jié)論。根據(jù)已有結(jié)論，得到SCAD回歸問題的穩(wěn)定點(diǎn)可保證良好的統(tǒng)計(jì)性質(zhì)。在分析的基礎(chǔ)上，我們設(shè)計(jì)了一種隨機(jī)坐標(biāo)下降方法(variable bregman stochastic coordinate descent, VBSCD)，這一方法能在內(nèi)存可控的前提下，很好地求解SCAD回歸問題，且該方法任意收斂點(diǎn)均是SCAD回歸模型的穩(wěn)定點(diǎn)。此外，在該算法每一次迭代的子問題中，更新變量數(shù)量可根據(jù)需求自由選擇，這一點(diǎn)很好地控制了算法每一次迭代中的內(nèi)存需求。本文通過計(jì)算實(shí)驗(yàn)說明這一方法在求解SCAD回歸問題時(shí)的有效性，決策人可根據(jù)隨機(jī)坐標(biāo)集合大小來控制子問題運(yùn)算所需要的內(nèi)存需求。

本文結(jié)構(gòu)如下：第1節(jié)對SCAD回歸問題的模型性質(zhì)進(jìn)行分析；第2節(jié)介紹求解非凸問題的隨機(jī)坐標(biāo)下降方法及其主要的理論結(jié)論，以及求解SCAD回歸的隨機(jī)坐標(biāo)下降算法流程；第3節(jié)進(jìn)行數(shù)值實(shí)驗(yàn)，說明本文方法的有效性；最后給出總結(jié)。

1 SCAD回歸問題模型的性質(zhì)分析

本節(jié)中考慮SCAD回歸模型，即

(7)

其中，稀疏懲罰項(xiàng)φλ(θ)為如下形式：

(8)

(9)

命題1：

‖▽f(x)-▽f(y)‖≤L‖x-y‖，?x,y∈domf,

(2)懲罰項(xiàng)g(β)是semi-convex函數(shù)，即存在常數(shù)ρ≥0，

證明：

(1)根據(jù)定義有

(2)這里可直接對φλ(θ)進(jìn)行分析，當(dāng)θ>0時(shí)，φλ(θ)的次微分為

因此有

學(xué)習(xí)情境的設(shè)定是根據(jù)實(shí)際工作過程提出的，這要求教師有工程意識、過程意識及全局意識，要求教師具有行業(yè)從業(yè)經(jīng)驗(yàn)。因此，制約地方新建本科院校向應(yīng)用技術(shù)大學(xué)轉(zhuǎn)型的瓶頸應(yīng)該是具有行業(yè)背景的師資。

此外，根據(jù)[23]的結(jié)論，我們可知SCAD回歸模型的任意穩(wěn)定點(diǎn)與全局最優(yōu)點(diǎn)具有基本一致的統(tǒng)計(jì)性質(zhì)。通過上述分析，不難看出，算法可求得SCAD回歸問題的穩(wěn)定點(diǎn)，即可得到原始參數(shù)的較好估計(jì)。下一節(jié)中將應(yīng)用本節(jié)中分析到的SCAD回歸模型的損失函數(shù)是導(dǎo)數(shù)Lipschitz、懲罰函數(shù)是semi-convex的分析結(jié)論。

2 求解SCAD回歸問題的隨機(jī)坐標(biāo)下降方法

2.1 求解非凸問題的隨機(jī)坐標(biāo)下降方法

最近Zhao和Zhu(2019)給出了求解如下非凸問題的隨機(jī)坐標(biāo)下降方法，針對的問題如下：

(8)

假設(shè)1：

(i)f是一個(gè)光滑函數(shù)，domf是凸的。f的梯度▽f是L-Lipschitz的；

(iii)F是水平有界的，即對任意r∈R，集合{x∈Rn|F(x)≤r}是有界的。

針對非凸問題(8)的隨機(jī)坐標(biāo)下降方法(variable bregman stochastic coordinate descent, VBSCD)流程如下：

選擇初始解x0∈Rn

fork=0,1,…, do

從{1,…,N}中等概率選擇i(k)

(9)

end for

在算法流程中D為Bregman函數(shù)，即D(x,y)=K(y)-[K(x)+〈▽K(x),y-x〉]，在Bregman函數(shù)中，K為核函數(shù)，為步長參數(shù)。核函數(shù)K和步長參數(shù)服從假設(shè)2。

假設(shè)2：

(i)K是m-強(qiáng)凸函數(shù)，且其梯度▽K是M-Lipschitz的。

在非凸問題(8)滿足假設(shè)1，算法VBSCD滿足假設(shè)2的情況下，得到了如下結(jié)論：算法VBSCD產(chǎn)生的任意收斂點(diǎn)是函數(shù)的一個(gè)穩(wěn)定點(diǎn)。

結(jié)合本節(jié)中算法性質(zhì)和第二節(jié)命題1中SCAD回歸模型性質(zhì)，利用Zhao和Zhu最近提出的VBSCD算法可以求解SCAD回歸問題。

2.2 求解SCAD回歸問題的隨機(jī)坐標(biāo)下降算法

求解SCAD回歸問題的隨機(jī)坐標(biāo)下降方法如下：

選擇初始解x0∈Rn

fork=0,1,…, do

(10)

從{I1,…,IN}中等概率選擇I(k)，對任意j∈I(k)，計(jì)算

(11)

end for

其中，子問題(11)有如下閉合表達(dá)式(見Fan和Li, 2001)：

(12)

其中，rk=βjk+

通過上述算法流程，不難看出本文給出的VBSCD算法，在算法每一次更新中僅隨機(jī)選擇一部分變量進(jìn)行更新，這一點(diǎn)很好地適應(yīng)了大數(shù)據(jù)時(shí)代求解SCAD回歸問題時(shí)數(shù)據(jù)規(guī)模大、分布在不同區(qū)域的需求。更新變量個(gè)數(shù)與所選擇隨機(jī)坐標(biāo)集合I(k)的大小緊密相關(guān)。每一次迭代中，所需要的內(nèi)存大小也和所選擇隨機(jī)坐標(biāo)集合I(k)的大小緊密相關(guān)。下一節(jié)，我們將通過仿真實(shí)驗(yàn)，說明每一次迭代中選擇隨機(jī)坐標(biāo)集合I(k)的大小與子問題求解內(nèi)存需求、算法總迭代次數(shù)之間的關(guān)系。

3 數(shù)值實(shí)驗(yàn)

本節(jié)對兩組仿真數(shù)據(jù)進(jìn)行數(shù)值實(shí)驗(yàn)，通過數(shù)值實(shí)驗(yàn)說明：

(1)本文隨機(jī)坐標(biāo)下降算法可有效求解SCAD回歸問題；

(2)本文隨機(jī)坐標(biāo)下降算法對不同規(guī)模算例求解時(shí)，所需迭代次數(shù)相對穩(wěn)定；

(3)本文隨機(jī)坐標(biāo)下降算法在不同變量分組情況下，算法迭代到穩(wěn)定點(diǎn)所需的迭代回合數(shù)(epoch)相對穩(wěn)定；

(4)本文隨機(jī)坐標(biāo)下降算法隨著變量分組的增加[隨機(jī)坐標(biāo)集合I(k)變小]，每一次迭代中子問題計(jì)算內(nèi)存需求減少。

3.1 數(shù)據(jù)集與計(jì)算環(huán)境

(1)仿真數(shù)據(jù)：樣本因素?cái)?shù)據(jù)X∈Rn×p的選擇服從N(0,1)高斯分布，構(gòu)造真實(shí)系數(shù)向量β*∈Rp，維數(shù)為p，稀疏度為s，β*中隨機(jī)的s個(gè)元素為非0元素，其余元素為0。非0元素的值服從N(0,1)高斯分布。反應(yīng)變量y=Xβ*+ε，其中各元素的噪聲向量獨(dú)立同分布，服從N(0,σ2)，用σ來控制噪聲大小。

根據(jù)上述規(guī)則，生成兩組仿真實(shí)驗(yàn)，分別為σ=0.01、n=100、p=1000、s=10；以及σ=0.01、n=200、p=2000、s=10。

(2)計(jì)算環(huán)境：Intel Core i5-6200U CPUs (2.40GHz)，8.00 GB RAM

3.2 實(shí)驗(yàn)結(jié)論

圖2和圖4中實(shí)線表示將變量分為100組時(shí)目標(biāo)函數(shù)的下降情況，點(diǎn)線表示將變量分為50組時(shí)目標(biāo)函數(shù)的下降情況，點(diǎn)劃線表示將變量分為10組時(shí)目標(biāo)函數(shù)的下降情況，虛線表示將變量分為5組時(shí)的目標(biāo)函數(shù)下降情況。

圖3和圖5中圓圈(o)表示原始系數(shù)，星號(*)表示通過SCAD回歸模型計(jì)算出的系數(shù)。

圖2 n=100、p=1000、s=10時(shí)變量劃分塊數(shù)與算法下降效果之間的關(guān)系

圖3 n=100、p=1000、s=10時(shí)信息恢復(fù)效果

圖8表示對n=100、p=1000、s=10問題將變量分為1、2,5、10、50、100組時(shí)，求解子問題所需的內(nèi)存大小。圖9表示對n=200、p=2000、s=10問題將變量分為1、2、5、10、50、100組時(shí)，求解子問題所需的內(nèi)存大小。

圖4 n=200、p=2000、s=10時(shí)變量劃分塊數(shù)與算法下降效果之間的關(guān)系

圖5 n=200、p=2000、s=10時(shí)信息恢復(fù)效果

圖6 n=200、p=2000、s=10時(shí)變量劃分塊數(shù)、算法迭代回合(epoch)、目標(biāo)函數(shù)下降效果之間的關(guān)系

圖7 n=200、p=2000、s=10時(shí)變量劃分塊數(shù)、算法迭代回合(epoch)、目標(biāo)函數(shù)下降效果之間的關(guān)系

圖8 n=100、p=1000、s=10算例，當(dāng)把變量分成1、2、5、10、50、100組時(shí)子問題內(nèi)存需求

圖9 n=200、p=2000、s=10算例，當(dāng)把變量分成1、2、5、10、50、100組時(shí)子問題內(nèi)存需求

通過對圖2至圖9的分析有如下結(jié)論：

(1)本文隨機(jī)坐標(biāo)下降算法可有效求解SCAD回歸問題：從圖2和圖4可以看出算法在迭代3000次左右時(shí)，即可收斂到較為穩(wěn)定的解。從圖3和圖5可以看出SCAD回歸具有很好的信息回歸效果。圖中圓圈(o)表示原始系數(shù)，星號(*)表示通過SCAD回歸模型計(jì)算出的系數(shù)。

(2)本文隨機(jī)坐標(biāo)下降算法對不同算例規(guī)模所需迭代次數(shù)相對穩(wěn)定：如圖2和圖4中，實(shí)線、點(diǎn)線、點(diǎn)劃線、虛線的對比可以得到，隨著拆分成的組(block)數(shù)量減少[隨機(jī)坐標(biāo)集合I(k)規(guī)模增加]，算法收斂到所需要結(jié)果的迭代次數(shù)減少。此外，對比圖2和圖4，可以得到不同規(guī)模算例，在分組數(shù)一致的情況下，迭代到所需要點(diǎn)的所需迭代次數(shù)基本相同。

(3)本文隨機(jī)坐標(biāo)下降算法在不同變量分組情況下，算法迭代到穩(wěn)定點(diǎn)所需的迭代回合數(shù)(epoch)相對穩(wěn)定：如圖6和圖7所示，在兩組算例中，目標(biāo)函數(shù)值隨迭代回合數(shù)增加而下降，其下降效率受變量分組數(shù)(隨機(jī)坐標(biāo)集合I(k)大小)影響不大。

(4)本文隨機(jī)坐標(biāo)下降算法隨著變量分組的增加[隨機(jī)坐標(biāo)集合I(k)變小]，每一次迭代中子問題計(jì)算內(nèi)存需求減少：如圖8和圖9中可以得到，隨著變量分組的增加[隨機(jī)坐標(biāo)集合I(k)變小]，算法子問題所需計(jì)算內(nèi)存減少。對比圖8和圖9，可以看出每組中變量的個(gè)數(shù)是子問題求解內(nèi)存需求的關(guān)鍵因素，在每組選擇變量個(gè)數(shù)相近的情況下，子問題的計(jì)算內(nèi)存需求相近。

4 總結(jié)

SCAD回歸問題，以其在處理高維數(shù)據(jù)中的近似無偏性，在大數(shù)據(jù)分析中得到廣泛應(yīng)用。但由于其懲罰項(xiàng)的非凸性，現(xiàn)有的隨機(jī)坐標(biāo)下降方法難以處理這一問題。本文首先對SCAD回歸模型進(jìn)行分析，得到SCAD回歸模型的損失函數(shù)是導(dǎo)數(shù)Lipschitz，懲罰函數(shù)是semi-convex的，并根據(jù)已有結(jié)論得出SCAD回歸模型的任意穩(wěn)定點(diǎn)有良好的統(tǒng)計(jì)性質(zhì)?；谶@些分析，本文介紹一種新的隨機(jī)坐標(biāo)下降方法，可應(yīng)用于求解大規(guī)模SCAD回歸問題，算法的任意收斂點(diǎn)均是SCAD回歸模型的穩(wěn)定點(diǎn)。這一方法能很好地求解SCAD回歸問題。本文通過計(jì)算實(shí)驗(yàn)說明這一方法在求解SCAD回歸問題時(shí)的有效性，同時(shí)說明算法在不同變量分組情況下，算法迭代到穩(wěn)定點(diǎn)所需的迭代回合數(shù)(epoch)相對穩(wěn)定。隨著變量分塊數(shù)的增加，單次迭代中計(jì)算的內(nèi)存需求減少。這一方法可被廣泛應(yīng)用于大數(shù)據(jù)分析和決策實(shí)踐當(dāng)中。

在大數(shù)據(jù)分析工作里，求解非凸非光滑的隨機(jī)坐標(biāo)下降算法的計(jì)算復(fù)雜度分析是個(gè)重要問題，Zhao和Zhu證明了在水平集次微分誤差界條件下，隨機(jī)坐標(biāo)下降算法可以達(dá)到線性收斂，而Li和Pong證明了SCAD回歸問題的目標(biāo)函數(shù)滿足該誤差界條件。因文章篇幅問題，本文算法的計(jì)算復(fù)雜度分析不在此處做進(jìn)一步討論。