石躍勇　鄧世容　焦雨領(lǐng)　徐德義

[摘 要]以統(tǒng)計(jì)學(xué)習(xí)（statistical learning）課程中基于高維線性模型的LASSO正則化為例，項(xiàng)目組給出了其基于原始對(duì)偶有效集（primal dual active set）方法的求解過程，并進(jìn)行了相應(yīng)的實(shí)證分析，期望培養(yǎng)學(xué)生對(duì)統(tǒng)計(jì)學(xué)習(xí)課程的學(xué)習(xí)興趣和應(yīng)用新近文獻(xiàn)方法進(jìn)行高維數(shù)據(jù)分析的實(shí)踐經(jīng)驗(yàn)。

[關(guān)鍵詞]原始對(duì)偶有效集方法； 高維稀疏； LASSO正則化；實(shí)證分析

[中圖分類號(hào)] G642 [文獻(xiàn)標(biāo)識(shí)碼] A [文章編號(hào)] 2095-3437（2019）05-0101-03

一、引言

大數(shù)據(jù)的興起和普遍，引起了學(xué)界和業(yè)界日益廣泛的重視[1-2]。在大數(shù)據(jù)時(shí)代，如在氣候研究、網(wǎng)絡(luò)交易、圖像處理等諸多領(lǐng)域，人們可以輕松獲得PB級(jí)、EB級(jí)、ZB級(jí)甚至更大量級(jí)的海量數(shù)據(jù)。直觀上看，大數(shù)據(jù)具有樣本量大或維度高等顯著特點(diǎn)。如何對(duì)大數(shù)據(jù)進(jìn)行建模和分析，提取蘊(yùn)藏其中的有用信息，進(jìn)而指導(dǎo)實(shí)際生產(chǎn)生活有著重要的現(xiàn)實(shí)意義。然而，大數(shù)據(jù)由于存在Volume（大量）、Velocity（高速）、Variety（多樣）、Value（價(jià)值）等特性[2]，給相應(yīng)的存儲(chǔ)、分析和應(yīng)用帶來了巨大的挑戰(zhàn)。特別地，有一類特殊的大數(shù)據(jù)，其具有很高的維度（可以達(dá)到數(shù)萬、數(shù)十萬甚至更高），而樣本量往往只有數(shù)十至數(shù)百，通常還具有復(fù)雜的結(jié)構(gòu)，我們稱之為高維數(shù)據(jù)。數(shù)據(jù)的高維度一方面可以反映研究對(duì)象更多的特征和屬性，可能使得相關(guān)的數(shù)據(jù)分析更能接近事物的本原，但另一方面也會(huì)產(chǎn)生“維數(shù)災(zāi)難”（ the curse of dimensionality）[3]，可能會(huì)帶來較為嚴(yán)重的過擬合（over-fitting）或不適定性（ill-posed）等諸多問題。所幸的是，在實(shí)際問題中，高維數(shù)據(jù)通常具有所謂的稀疏性結(jié)構(gòu)，即影響某種結(jié)果的因素有很多，但關(guān)鍵影響因素卻較少。基于稀疏性假設(shè)，人們可以借助于統(tǒng)計(jì)學(xué)習(xí)（機(jī)器學(xué)習(xí)）中的正則化方法[3]對(duì)高維數(shù)據(jù)進(jìn)行建模、分析和應(yīng)用。

正則化方法是處理高維數(shù)據(jù)的利器，也是統(tǒng)計(jì)學(xué)習(xí)（機(jī)器學(xué)習(xí)）等學(xué)科的重要內(nèi)容，其模型、理論和計(jì)算等方面的研究是人們關(guān)注的焦點(diǎn)[1，3]。由于實(shí)際問題中高維數(shù)據(jù)不斷地、大量地涌現(xiàn)，當(dāng)下一個(gè)重要的發(fā)展趨勢是：人們?cè)谶M(jìn)行數(shù)據(jù)分析和統(tǒng)計(jì)計(jì)算時(shí)，越來越注重快速算法的導(dǎo)向，即兼顧統(tǒng)計(jì)性質(zhì)（如無偏性、相合性等[4]）和數(shù)學(xué)收斂的同時(shí)，在保證精度（估計(jì)誤差、預(yù)測誤差等）相當(dāng)?shù)那疤嵯?，希望設(shè)計(jì)高效的算法使得計(jì)算速度越快越好。這一趨勢在大數(shù)據(jù)高速發(fā)展的大背景下日益顯著和越發(fā)重要，可在此基礎(chǔ)上進(jìn)一步研究適用于分布式環(huán)境下的稀疏化學(xué)習(xí)理論和算法。

統(tǒng)計(jì)學(xué)習(xí)是一套以理解數(shù)據(jù)為目的的龐大工具集[5]，其研究對(duì)象和內(nèi)容與模式識(shí)別、機(jī)器學(xué)習(xí)、數(shù)據(jù)挖掘等學(xué)科有大量的交叉和重疊，是數(shù)據(jù)科學(xué)的重要組成部分，廣泛應(yīng)用于社會(huì)生活的各個(gè)方面。在當(dāng)前大數(shù)據(jù)研究蓬勃發(fā)展的大背景下，統(tǒng)計(jì)學(xué)習(xí)的課堂教學(xué)和學(xué)科發(fā)展的研究顯得尤為重要。鑒于統(tǒng)計(jì)學(xué)習(xí)學(xué)科具有日新月異、發(fā)展迅速等特點(diǎn)，在統(tǒng)計(jì)學(xué)習(xí)的教學(xué)過程中，一個(gè)重要的趨勢是，除了講授經(jīng)典的理論和方法，還應(yīng)當(dāng)引導(dǎo)學(xué)生閱讀新近的文獻(xiàn)，鼓勵(lì)學(xué)生運(yùn)用文獻(xiàn)中新的建模和計(jì)算方法開展數(shù)據(jù)分析與應(yīng)用，并結(jié)合Julia、Matlab、Python或R等軟件進(jìn)行相應(yīng)的數(shù)值模擬和實(shí)證分析。特別地，R （https：//www.r-project.org/）作為一款自由的編程軟件，提供大量最新的開源的方法和算法程序包，社區(qū)成熟，參與者眾多，在Linux、Mac、Windows平臺(tái)下均可方便地進(jìn)行安裝和使用，深受學(xué)界歡迎。學(xué)生在統(tǒng)計(jì)學(xué)習(xí)課程中學(xué)習(xí)使用 R 軟件，熟練掌握R技能，并將所學(xué)方法用R編程實(shí)現(xiàn)，有助于其更為深刻地聯(lián)系理論和理解問題，且能夠更為有效地拓展實(shí)際應(yīng)用，提高學(xué)習(xí)的興趣和實(shí)戰(zhàn)的能力[6]。

線性模型是研究其他統(tǒng)計(jì)模型的基礎(chǔ)，線性模型正則化[3，5]的理論和計(jì)算是統(tǒng)計(jì)學(xué)習(xí)的基本內(nèi)容和研究熱點(diǎn)，其中又以LASSO正則化（即[?1]正則化）模型尤為重要，因?yàn)長ASSO正則化模型的求解屬于凸優(yōu)化問題，凸性可以保證模型的解具有良好的理論分析性質(zhì)和數(shù)值收斂結(jié)果。LASSO正則化在壓縮感知（信號(hào)、圖像處理）和變量選擇（高維特征降維）的建模和計(jì)算中均得到了廣泛的應(yīng)用。本文使用文獻(xiàn)[7]中發(fā)展的原始對(duì)偶有效集方法，對(duì)高維線性模型的LASSO正則化問題進(jìn)行了應(yīng)用，對(duì)一組高維基因微陣列數(shù)據(jù)進(jìn)行了變量選擇實(shí)證分析，并將計(jì)算結(jié)果與著名的R統(tǒng)計(jì)軟件包glmnet進(jìn)行比較，在開闊學(xué)生視野、拓寬研究思路、豐富學(xué)習(xí)內(nèi)容的同時(shí)，期望加深學(xué)生對(duì)正則化模型及其計(jì)算的理解。

二、模型與方法

其中[y∈Rn]是響應(yīng)變量，[X∈Rn×p]是設(shè)計(jì)矩陣，[β∈Rp]是待估系數(shù)向量，[ò∈Rn]是噪聲向量。不失一般性，我們考慮噪聲滿足正態(tài)（高斯）分布[ò～N（0， σ2In）]，這里[In]表示[n×n]單位矩陣。進(jìn)一步地，我們?cè)O(shè)定數(shù)據(jù)滿足高維情形[p>n]，并假設(shè)向量[β]的非零分量個(gè)數(shù)是稀疏的。由于[p>n]，傳統(tǒng)的統(tǒng)計(jì)方法（如最小二乘）對(duì)求解模型（1）的參數(shù)估計(jì)已不再適用，因?yàn)榈玫降慕獠皇俏ㄒ坏??；谙∈栊约僭O(shè)，我們考慮高維線性模型的LASSO正則化如下：

其中[?]和[‖?‖1]分別表示歐式范數(shù)（[?2]范數(shù)）和[?1]范數(shù)，[λ>0]為調(diào)節(jié)參數(shù)（正則化參數(shù)）。問題的求解可以直接調(diào)用R軟件中廣受歡迎的glmnet軟件包[3]，該包基于坐標(biāo)下降（coordinate descent）算法，可以求解線性模型和logistic模型的[?1]正則化問題，具體使用方法可以查閱軟件包的參考手冊(cè)或教材[5]第6章第6.6節(jié)實(shí)驗(yàn)2的內(nèi)容。

文獻(xiàn)[7]發(fā)展了一種原始對(duì)偶有效集方法并結(jié)合一種連續(xù)化策略求解問題（2），方法英文簡記為PDASC。該方法在壓縮感知（一維信號(hào)和二維磁共振圖像）中進(jìn)行了有效的應(yīng)用，其主要步驟概括為：

1.基于問題中目標(biāo)泛函[Qβ]的坐標(biāo)極小值[β*]（原始變量）給出對(duì)偶變量[d*=XT（y-Xβ*）]，且對(duì)LASSO模型給出關(guān)系表達(dá)式[β*=Sλ（β*+d*）]，其中[Sλ（t）=max（|t|-λ，0）sgn（t）]為軟閾值算子；

2.對(duì)LASSO模型有[|β*j+d*j|>T*?β*j≠0]成立，其中[T*=λ]，從而可以構(gòu)造有效集[A*=j：|β*j+d*j|> T*]和非有效集[I*=（A*）c]；

3.設(shè)[βk]和[dk]分別逼近[β*]和[d*]，基于[βk]和[dk]得到[Ak+1=j：|βkj+dkj|> T*]和[Ik+1=（Ak+1）c]，根據(jù)構(gòu)造知[Ak+1]和[Ik+1]分別逼近[A*]和[I*]，再基于[Ak+1]和[Ik+1]求解一個(gè)低維最小二乘問題得到[βk+1]和[dk+1]，對(duì)[k]作循環(huán)，設(shè)置停機(jī)準(zhǔn)則（[Ak=Ak+1]或[k≥K]，[K]為某個(gè)正整數(shù)），得到一個(gè)給定[λ]的估計(jì)[β（λ）=βk+1]；

4.考慮一串[λ]格子點(diǎn)序列[λss]，其中[λs=λ0ρs]，[λ0=||XTy||∞]，[0<ρ<1]，對(duì)于每個(gè)[λs]用步驟3得到估計(jì)[β（λs）]，再結(jié)合連續(xù)化策略（continuation strategy）得到整個(gè)解的路徑[{β（λs）}s]，最后根據(jù)BIC或HBIC等調(diào)節(jié)參數(shù)選擇準(zhǔn)則選擇最優(yōu)的[λ]和最優(yōu)的解[β（λ）]。

可以看出，PDASC方法操作簡便、實(shí)施簡單，關(guān)鍵之處是同時(shí)利用原始變量和對(duì)偶變量的信息確定一個(gè)有效集（相比之下，經(jīng)典的LARS方法僅僅使用了對(duì)偶變量的信息），然后在有效集上求解一個(gè)低維問題。PDAS算法具有局部超線性的收斂速度，而連續(xù)化策略可以使得問題的求解具有良好的初值（因?yàn)镻DAS方法屬于Newton型算法，其收斂性依賴于初值的選?。?，再配合PDAS算法的局部超線性收斂性，可以使得整體的PDASC步驟計(jì)算速度非?？?。特別地，在理論上，文獻(xiàn)[7]在LASSO正則化下證明了PDAS的局部一步收斂性（推廣了現(xiàn)有的局部超線性收斂結(jié)果），并且在壓縮感知框架下證明了PDASC對(duì)于LASSO模型的全局收斂性。更多關(guān)于PDASC方法的理論細(xì)節(jié)和實(shí)施步驟（算法偽代碼、調(diào)節(jié)參數(shù)選擇準(zhǔn)則等）參見文獻(xiàn)[7]。在高通量生物醫(yī)學(xué)研究中，例如基因關(guān)聯(lián)研究和基因表達(dá)研究，基于正則化的變量選擇方法被大量應(yīng)用，而設(shè)計(jì)具有較高精度的快速算法求解相應(yīng)正則化問題是應(yīng)用的關(guān)鍵。在下一節(jié)，我們運(yùn)用PDASC方法對(duì)DNA微陣列基因表達(dá)數(shù)據(jù)在高維LASSO模型框架下展開計(jì)算和分析，以期望獲得較為合理的具有一定生物學(xué)意義的結(jié)果。

三、實(shí)證分析

根據(jù)文獻(xiàn)[8]，我們對(duì)ISLR軟件包[5]中的NCI60微陣列數(shù)據(jù)進(jìn)行實(shí)證分析，主要目標(biāo)是在高維LASSO模型框架下比較glmnet和PDASC兩種方法的計(jì)算速度和精度。NCI60數(shù)據(jù)包含了來自64個(gè)癌細(xì)胞系的6830個(gè)基因的表達(dá)水平，關(guān)于該數(shù)據(jù)的更多信息可在網(wǎng)站http：//genome-www.stanford.edu/nci60/ 上獲取。我們的研究目的是考察第1個(gè)基因和余下所有基因的關(guān)系。在線性模型（1）下，該數(shù)據(jù)的響應(yīng)變量[y]是長度為64的數(shù)值向量，其給出了第1個(gè)基因的表達(dá)水平；設(shè)計(jì)矩陣[X]是維度為[ 64 × 6829]的矩陣，其表示了剩余的6829個(gè)基因的表達(dá)水平。故該數(shù)據(jù)有[n=64]，[p=6829]，滿足[p>n]，從而是一個(gè)高維數(shù)據(jù)，于是我們?cè)谡齽t化模型（2）下求解問題。對(duì)于glmnet和PDASC兩種計(jì)算方法，我們均采用文獻(xiàn)[9]中的voting準(zhǔn)則選擇調(diào)節(jié)參數(shù)[λ]。對(duì)于兩種方法（Method），我們采用的比較指標(biāo)包括計(jì)算時(shí)間（Time，單位：秒）、模型大?。∕S，即非零元素個(gè)數(shù)）、預(yù)測誤差（PE）和選擇的基因（Gene，即設(shè)計(jì)矩陣[X]對(duì)應(yīng)的列），其中Time為速度指標(biāo)，其余為精度指標(biāo)，結(jié)果如表1所示。由表1可知，PDASC的PE比glmnet大，但計(jì)算速度比glmnet快，且選擇的基因個(gè)數(shù)比glmnet少。由表1還可知，PDASC選擇的基因集合是glmnet選擇基因集的真子集，這說明了PDASC基因集具有一定的顯著性和代表性。我們?cè)诒?中列出了兩種方法相同基因的估計(jì)系數(shù)，發(fā)現(xiàn)雖然兩種方法給出的估計(jì)系數(shù)絕對(duì)值不同，但符號(hào)相同，說明二者表達(dá)的生物學(xué)意義是一樣的。特別地，表2中PDASC給出的Gene估計(jì)系數(shù)絕對(duì)值要比glmnet給出的大，注意到PDASC基因集比glmnet基因集小，從而在一定意義上說明PDASC方法可以更為集中而精確地選擇到顯著的自變量集。 值得注意的是，為了使結(jié)果更具有說服力，還應(yīng)該使用兩種方法對(duì)NCI60數(shù)據(jù)分析作交叉驗(yàn)證（cross validation），以獲得更為穩(wěn)健的速度和精度結(jié)果，可參考教材[5]第6章第6.6節(jié)實(shí)驗(yàn)2的相關(guān)內(nèi)容。

四、結(jié)語

本文對(duì)統(tǒng)計(jì)學(xué)習(xí)課程中基于高維線性模型的LASSO正則化的計(jì)算問題進(jìn)行了拓展訓(xùn)練，結(jié)合統(tǒng)計(jì)學(xué)習(xí)學(xué)科發(fā)展迅速的特點(diǎn)，通過引入文獻(xiàn)中一種新的原始對(duì)偶有效集方法對(duì)高維基因微陣列數(shù)據(jù)作實(shí)證分析，并與著名的glmnet包的計(jì)算結(jié)果進(jìn)行了比較。

通過本文內(nèi)容的介紹和實(shí)證，我們希望拋磚引玉，引導(dǎo)同學(xué)們?cè)趯W(xué)習(xí)教材內(nèi)容的同時(shí)，更加注重將課本知識(shí)與最新的科研方法相結(jié)合，多檢索、多思考、多動(dòng)手，為以后在學(xué)界從事科研教學(xué)或在業(yè)界進(jìn)行數(shù)據(jù)挖掘工作打下良好的基礎(chǔ)。

與此同時(shí)，除了學(xué)生們的“學(xué)”，我們認(rèn)為大數(shù)據(jù)時(shí)代下統(tǒng)計(jì)學(xué)習(xí)學(xué)科的快速發(fā)展也對(duì)老師們的“教”提出了新的更高的要求。廣大教師必須要更加切實(shí)地轉(zhuǎn)變教學(xué)思想，目標(biāo)和眼界要聚焦于領(lǐng)袖學(xué)校和行業(yè)巨頭的發(fā)展動(dòng)態(tài)，不斷適應(yīng)社會(huì)發(fā)展，及時(shí)更新統(tǒng)計(jì)學(xué)習(xí)教材，運(yùn)用多種教學(xué)方法和軟件手段來培養(yǎng)學(xué)生的學(xué)習(xí)興趣和操作經(jīng)驗(yàn)，鼓勵(lì)學(xué)生更加注重案例分析和社會(huì)實(shí)踐，鞏固學(xué)生掌握傳統(tǒng)經(jīng)典知識(shí)的同時(shí)，把握好“課堂所學(xué)”與“社會(huì)所需”之間的緊密聯(lián)系，激發(fā)學(xué)生大膽探索適應(yīng)大數(shù)據(jù)時(shí)代發(fā)展的新思想、新理論、新方法和新應(yīng)用。

[ 參 考 文 獻(xiàn) ]

[1] 焦雨領(lǐng).稀疏約束下反問題理論與算法的研究[D].武漢：武漢大學(xué)博士學(xué)位論文，2014.

[2] 徐德義，林志恒.對(duì)大數(shù)據(jù)時(shí)代大學(xué)統(tǒng)計(jì)教學(xué)的認(rèn)識(shí)與思考[J].大學(xué)教育，2015（11）：183-184.

[3] Hastie， T， Tibshirani， R， Friedman， J. Elements of statistical learning： data mining， inference， and prediction （second edition）[M]. Springer， 2009.

[4] 鄭明，陳子毅，汪嘉岡.數(shù)理統(tǒng)計(jì)講義[M].上海：復(fù)旦大學(xué)出版社，2006.

[5] [美]詹姆斯（James， G.）等著，王星等譯. 統(tǒng)計(jì)學(xué)習(xí)導(dǎo)論——基于R應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2015.

[6] 鄧世容，石躍勇.蒙特卡洛方法在概率論與數(shù)理統(tǒng)計(jì)教學(xué)中的應(yīng)用[J].科教導(dǎo)刊（上旬刊），2018（1）：111-112.

[7] Fan Q， Jiao Y， Lu X. A primal dual active set algorithm with continuation for compressed sensing[J]. IEEE Transactions on Signal Processing， 2014， 62（23）：6276-6285.

[8] Shi Y， Cao Y， Yu J， Jiao Y. High-dimensional variable selection with the generalized SELO penalty [J].Jounal of Mathematics， 2018， 38（6）， 900-998.

[9] Huang J， Jiao Y， Lu X， Zhu L. Robust decoding from 1-bit compressive sampling with ordinary and regularized least squares [J]. SIAM Journal on Scientific Computing， 2018， 40（4）：A2062-A2086.

[責(zé)任編輯：林志恒]