統計回歸模型及其優化算法綜述

方升，梁飛豹，劉勇進

(福州大學數學與統計學院，福建 福州 350108)

0 引言

由于正則化回歸模型可以用于變量選擇、 參數估計和預測，因此，其為當今處理高維數據的一種重要手段. 在過去的數十年中，高維回歸分析有著持續可觀的進展，具有代表性的統計模型有Lasso[1]、 Dantzig選擇器[2]，以及具有優異統計特性的凹正則項如SCAD[3]、 MCP[4]和截斷l1[5]等. 其中，因平方損失函數優良的數學特性，正則化最小二乘回歸模型得到了廣泛的應用. 但是，正則化最小二乘回歸模型對于隨機誤差服從非正態分布的情況并不具有魯棒性. 魯棒回歸模型，如分位數回歸模型(如最小絕對離差回歸)，可以用于處理誤差服從重尾分布的高維數據的情況. 歷年來，學者在如何提高正則化回歸模型的參數估計的準確性、 確保回歸模型的魯棒性、 降低變量估計的錯誤發現率和提高模型問題的計算效率等方面都有顯著的進展. 更多關于正則化回歸模型的綜述可參考文獻[6-7]. 針對大規模數據下的正則化回歸模型問題的優化算法[8-9]和算法軟件包也不斷涌現，相應的優化理論也得到不斷完善和發展，無論是計算簡便的一階算法[10]，如交替方向乘子法，還是近年來不斷崛起的二階算法，如半光滑牛頓增廣拉格朗日方法等，在求解大規模優化問題上都有良好的數值表現. 本研究主要回顧一些經典的正則化回歸模型，并對歷年來的優化算法進行了簡要的綜述.

1 正則項與優化算法

稀疏估計是當今處理高維數據的重要技術，一般通過求解正則化回歸模型問題進行變量選擇和稀疏向量的估計. 正則化是對特定的損失函數加上一個正則項，該正則項是關于模型復雜度的懲罰函數，從而使得問題的解具有稀疏性.

1.1 凸正則項

l2范數正則項最早由Tikhonov[11]提出，而后由Hoerl等[12- 13]將l2-正則化回歸模型引入統計學，并稱之為嶺回歸模型：

該模型用l2范數代替最優子集選擇中的l0范數，是最早的凸正則模型.

Candes等[2]于2007年提出了Dantzig選擇模型：

該模型可進一步轉化為線性規劃問題. 該估計量滿足非漸近Oracle不等式，具有正交變換不變性等許多優良統計特性.

稀疏度的概念和l1范數正則項首先出現在Donoho等[16]的工作中，Chen[17]使用l1-正則回歸模型，即基追蹤： 從過飽和字典中提取基. 受到非負絞刑估計等的啟發，Tibshirani[1]提出了最小絕對收縮與選擇算子模型(Lasso)進行變量選擇：

因Lasso模型便于求解，該模型得到了廣泛的應用和關注. 關于β的Lasso估計取決于預先知道噪聲的標準偏差，然而當特征維數很大特別是在特征維數遠大于樣本個數時，估計噪聲的標準偏差并不容易. 為克服這個缺點，Belloni等[18]提出了方根Lasso (srLasso) 模型來估計β：

該模型不需要知道或提前估計噪聲的標準偏差. 在一定條件下，即使噪聲的標準偏差未知，srLasso估計量仍具有極小極大最優收斂性[19- 20]. 由Sun等[21]提出的scaled Lasso估計與srLasso估計基本是等價的，而且，Xu等[22]證明了srLasso模型等價于在一定約束條件下的魯棒線性回歸問題.

2006年，Zou等[23]對Lasso進行改進并提出了自適應Lasso模型 (adaptive Lasso)，該模型具有Oracle性質，其具體形式如下：

其中:w∈Rp是已知的權重向量.

當特征(協變量)個數p大于樣本個數n時，Lasso至多只能選擇n個協變量(即使可能存在更多與因變量有關的協變量)； 當協變量間兩兩相關性很高時，Lasso只能選出其中一個協變量，無法選擇出其他相關協變量，然而在實際應用中，更傾向于找出所有相關變量 (如： 找出與某種疾病有關的所有變量)，而且只找出一個變量將會提高預測誤差，使得模型魯棒性更差； 另外，當n>p且協變量間有強相關性時，嶺回歸模型在變量選擇方面的表現要優于Lasso. 基于Lasso的上述缺陷，2005年，Zou等[24]結合Lasso和嶺回歸模型提出了彈性網絡(elastic net)正則化模型：

在實際應用中，有時不同協變量之間具有相關性，此時傾向于對協變量進行分組，進而研究因變量與不同組協變量之間的關系. 比如，在基因表達中，一個通路中的基因是強相關的，與某項基因表達特性有關的一般是某個通路而不是某個單獨的基因. 2006年，Yuan等[25]提出了如下對整組變量進行選擇或舍棄的組Lasso模型：

組Lasso模型不能對組內個別變量進行選擇，然而實際情況中需要對組內的個別變量進行選擇. 如在基因表達學習中，癌癥可能與某條通路有關，但并不是這條通路中的所有基因都是起作用的，此時不僅需要挑選某組基因，還需挑選該組中的個別基因. Friedman等[26]提出了如下稀疏組Lasso模型：

其中:λ1,λ2≥0.當λ1=0時，上式為組Lasso模型； 當λ2=0時，上式為Lasso模型. 稀疏組Lasso模型不僅可以獲得組間的稀疏性，也能獲得組內變量的稀疏性. 另外，針對不同組中存在重疊變量的情況，學者提出了重疊組Lasso模型[27].

與稀疏組Lasso相似，Exclusive Lasso正則項由Zhou等[28]提出用于促進不同組內的稀疏度，并將該正則項應用于多任務學習，其具體形式如下:

Bondell等[29]提出一有效的對特征進行分組的稀疏模型-OSCAR模型：

文獻[30]將上述的OSCAR模型轉化為如下的SLOPE模型：

其中: |β|(i)為|β|的第i個最大元素, |β|(1)≥|β|(2)≥…≥|β|(p)，λ1≥λ2≥…≥λp≥0,λ1>0.關于SLOPE模型的統計特性可參考文獻[31].

Lasso模型針對對象是每一個變量個體，無法處理連續變量數據(如時間序列數據和圖像數據). 2005年，Tibshirani等[32]對Lasso模型進行擴展，增加了關于相鄰變量間的差異的正則項，提出了融合Lasso模型：

其中:λ1,λ2≥0,Bβ=[β1-β2,β2-β3, …,βp-1-βp]T.該模型不僅可以獲得稀疏解，保證數據的局部連續性，而且更容易檢測數據的內在結構并建立模型捕捉它.作為融合Lasso模型的擴展，clustered Lasso模型[33-34]具有如下形式:

1.2 非凸正則項

最優子集選擇最早由Akaike等[35-37]提出，其形式可看作如下l0-正則最小二乘回歸問題：

其中:l0范數是向量中的非零元素的個數，其為非凸非連續函數.

1993年橋回歸模型[38- 39]被提出，其具有下列形式：

其中: 0<γ.當0<γ<1時，該正則項又稱擬范數，此時正則項是非凸的.由于γ>1時無法產生稀疏解，因此一般考慮0<γ<1.

當線性回歸方程中的設計矩陣X滿足一定的條件如受限等距條件[40]、 強不可表示條件[41- 42]以及回歸系數具有稀疏性時，Lasso估計才具有較好的變量估計和選擇特性. 然而，即使在上述條件下，Lasso估計仍存在偏差. 2001年，Fan等[3]提出一個好的正則項產生的估計量需要具備以下三條性質： 無偏性、 稀疏性和連續性，但lq范數正則化最小二乘回歸所產生的估計量中沒有一個同時滿足這三個性質，于是他提出如下同時具備上述性質的正則項——SCAD函數：

其中:pλ(βi)具有下列形式：

其中:λ≥0和a>2. 較弱的條件下，SCAD正則項產生的估計量具有Oracle性質.

2010年，Zhang[5]提出了截斷l1正則項 (Capped-l1)：

該正則項對原點附近的值實行l1范數懲罰，對較大點的值用常數懲罰. 該正則項產生的估計量理論上具有無偏性和Oracle性質，可看作l0范數正則項的近似，也可看作簡單的SCAD正則項，但結果不如SCAD光滑.

2011年，Zhang[4]提出了極小極大凹正則項 (MCP):

其中:a>1. 理論上該估計是近乎無偏估計且具有Oracle性質.

1.3 優化算法

經典的求解無約束優化問題的算法有： 梯度法、 牛頓法、 共軛梯度法和信賴域法等. 梯度法，顧名思義是利用目標函數的一階信息—梯度進行迭代的優化算法，是最早、 最簡單的方法. 梯度法的計算形式雖然簡單，但是收斂速度只能達到線性，因此學者對梯度法進行加速可使其具有更快收斂速度[43]. 牛頓法利用目標函數的二階信息-Hessian矩陣構造下降方向，雖可達到局部二次收斂，但計算代價較大. 求解約束優化問題的優化算法主要有： 可行點方法、 交替方向乘子法和增廣拉格朗日方法等. 實際應用中，可將上述不同算法相結合從而更高效求解目標問題. 接下來，根據算法中是否涉及到目標函數的一階信息和二階信息對相關算法進一步分類詳細介紹.

1.3.1 一階算法

一階算法雖無法獲得高精度解，但因其計算簡便快速，且多數情況下并不是必需要獲得高精度解，所以一階算法廣泛應用于大規模優化問題的求解. 下面主要介紹兩種經典且重要的一階算法： 加速鄰近梯度法和交替方向乘子法. 更多其他一階算法的介紹，可參考文獻[10].

1) 加速鄰近梯度法. 加速鄰近梯度法(APG)由Nesterov[43]提出，之后Beck等[44]將該方法應用于l1范數正則最小化問題，并稱方法為快速軟閾值迭代算法(FISTA).

考慮如下目標問題：

其中, 函數g是正常的閉凸函數(可能是非光滑的)，函數f是連續可微的凸函數且是Lf-光滑即：

任意L>0，考慮F(x)在給定點y處的如下二次近似：

下面給出固定步長L的FISTA的算法框架，關于步長選取的其他方法詳見文獻[44].

算法1: 固定步長L的FISTA初始化: 選取L>0, 設置y0=x0, t0=1. 對k=0, 1, …迭代直至滿足停止準則S1計算: xk+1=Pro xgLyk-▽f(yk)L S2計算: tk+1=1+1+4t2k2S3計算: yk+1=xk+1+(tk-1)xk+1-xktk+1

Liu等[46]開發了軟件包SLEP，該軟件包運用加速鄰近梯度法求解一系列稀疏優化問題如Lasso、 融合Lasso、 組Lasso、 稀疏組Lasso、 稀疏邏輯回歸問題等，該軟件包得到了廣泛的應用.

下面給出ADMM算法求解該目標問題的具體框架.

算法2: 交替方向乘子法(ADMM)初始化: 給定ρ>0, 選取(z0, y0), 對k=0, 1, …迭代直至滿足停止準則S1計算: xk+1:=arg minx Lρ(x, zk; yk)S2計算: zk+1:=arg minz Lρ(xk+1, z; yk)S3計算: yk+1:=yk+ρ(Axk+1+Bzk+1-c)

ADMM算法計算簡便快速，廣泛應用于機器學習和統計學大規模優化問題的求解, 如Dantzig選擇模型問題[56]、 壓縮感知方面的問題[57]、 分位數回歸模型問題[58]、 支持向量機[59]等. 考慮到ADMM算法是一階算法，無法收斂于高精度解，常常將其作為算法的第一階段(熱啟動warm-start)用于獲得低精度的解從而加速后續算法的收斂率，可參考文獻[60-62].

1.3.2 二階算法

由于利用了問題的二階信息，二階算法的收斂速度優于一階算法并且可得到高精度解，但利用二階信息使得算法的計算代價增加，因此一般的二階算法往往只在尋求高精度解的時候才能超越一階算法. 下面介紹兩種高效的二階算法： 鄰近牛頓法和半光滑牛頓增廣拉格朗日方法.

① 鄰近牛頓法. 考慮如下目標問題是復合函數的形式：

(1)

其中:f是連續可微的凸函數;g是凸函數但可能是不可微的.

許多求解上述形式的目標問題的算法都可看作鄰近牛頓法(或又稱作連續二次近似方法SQA)，如投影牛頓型方法[9]、 鄰近擬牛頓方法[63- 64]、 坐標梯度下降法[65]等.

(2)

其中:Hk是正定矩陣，用以近似Hessian陣▽2f(xk).下面進一步給出鄰近牛頓法的具體算法框架：

算法3: 鄰近牛頓法初始化: 選取x0, 對k=0, 1, …迭代, 直至滿足停止條件 S1選取近似Hessian陣▽2f(xk)的正定矩陣Hk S2求解子問題(2)獲得方向dk=x︵k+1-xk即得到: dk←arg mind ▽f(xk)Td+12dTHkd+g(xk+d) S3 運用回溯線搜索得到步長tk S4更新xk+1←xk+tkdk

由于鄰近牛頓法的收斂性質與內迭代解的精確度相關，Lee等[71]對鄰近牛頓法中的內迭代算法設置了自適應停止條件并證明了鄰近牛頓法的收斂性. Byrd等[72]提出了非精確求解子問題的兩條停止準則從而保證了算法的全局收斂性并可通過相應參數對鄰近牛頓法的局部收斂率進行調節，內迭代采用的是基于象限方法(OBM:orthant-based-method). 2019年，Yue等[73]提出新的非精確鄰近牛頓法-IRPN，相較于之前的鄰近牛頓法，IRPN更靈活且具有更滿意的收斂性質. IRPN不需要一個精確的內迭代算法，且式(1)中的函數f不需要是強凸的，IRPN產生的解全局收斂于最優解，并且當目標問題具有Luo-Tseng EB性質時，根據算法參數的選取，局部收斂率是線性、 超線性或甚至可達二次. 由于IRPN算法的局部二次收斂率需要參數滿足與誤差上界值有關的相應條件，在實際過程中是較難估計的. Mordukhovich等[74]提出了全局收斂、 局部超線性收斂(甚至二次收斂)的鄰近牛頓法，該方法對問題的強凸性不作要求，在收斂性的理論保證和數值表現上都優于IRPN算法.

② 半光滑牛頓增廣拉格朗日方法. 半光滑牛頓增廣拉格朗日方法(SSNAL算法)無論尋找高精度解還是低精度解的數值表現都遠優于一階算法. 文獻[60-61, 75]受到求解半定規劃的兩層迭代算法框架的啟發提出了SSNAL算法，該算法的外層是鄰近點算法[76](PPA)或增廣拉格朗日方法[77](ALM)(運用增廣拉格朗日方法求解目標問題相當于對目標問題的對偶問題采用鄰近點方法求解)，內迭代算法是半光滑牛頓方法. 在較弱的假設條件下，鄰近點算法具有優異的收斂性質： 全局收斂和局部漸近超線性收斂[76-78]. 而且，內迭代算法—半光滑牛頓法[79]的收斂率是超線性甚至可達二次[75]. 此外，SSNAL算法進一步對目標問題內蘊的非光滑二階信息進行深入研究，通過充分利用問題的二階稀疏性(廣義Hessian陣為0-1對角矩陣)極大減少了內迭代的工作量使得整體算法的工作量與一階算法相當或甚至更少.

由于SSNAL算法優異的收斂性質和低計算代價，SSNAL算法可以高效獲得大規模問題的高精度解. 該算法的高效性和魯棒性在求解許多大規模優化問題中都得到了驗證，如Lasso[80]、 組圖Lasso[81]、 融合Lasso[82]、 聚集(clustered)Lasso[83]、 稀疏組Lasso[84]、 Oscar和Slope[30]、 支持向量機[85]、 exclusive lasso[86]、 大規模的線性規劃問題[62]、 半定規劃領域的問題[60- 61, 75]、 帶有非凸正則項的方根損失回歸問題[87]等.

接下來介紹SSNAL算法的具體細節，考慮如下目標問題：

相應地，其對偶問題形式如下：

其中:h*(·)和p*(·)分別是h(·)和p(·)的共軛函數.

假設h是連續可微的凸函數并且其梯度是Lipschitz連續；h是本質局部強凸函數.該假設條件在復合凸優化問題中是常見的，如線性回歸和邏輯回歸中的損失函數均滿足上述假設條件.并且，當損失函數不滿足上述條件時，仍有其他辦法可以使得SSNAL算法用于求解目標問題(P)，詳見文獻[86- 87].

對偶問題的拉格朗日函數如下：

l(y,z;x)=h*(y)+p*(z)-〈x,A*y+z-c〉

并且可得相應的增廣拉格朗日函數：

下面給出求解問題(D)的SSNAL算法的具體框架：

算法4: 求解問題(D)的非精確增廣拉格朗日方法初始化: 給定σ0>0, 選取(y0, z0, x0), k:=0S1計算(yk+1, zk+1)≈arg min{Ψk(y, z):=Lσk(y, z; xk)}S2計算 xk+1=xk-σk(A*yk+1+zk+1-c)S3更新 σk+1↑σ∞≤∞, k=k+1

算法4中S1非精確求解子問題的停止準則，可詳見文獻[80].

其中: Proxp(x)是鄰近映射, 其定義如下：

ψ(·)是強凸的且連續可微，其導數如下：

▽ψ(y)=0

(3)

由假設條件及文獻[88]中的引理12.60可得：h*是可微且梯度是局部Lipschitz連續. 定義如下多值映射:

定義：

V:=H+σAUA*

在SSNAL算法中，采用半光滑牛頓法求解方程(3)，具體算法框架如下：

算法5: 求解問題(3)的半光滑牛頓法初始化: 給定μ∈0, 12 , η∈(0, 1), τ∈(0, 1], δ∈(0, 1), 選取 y0∈int(dom h*), 對j=0, 1, …, 迭代: S1選取Vj∈?︵2ψ(yj), 精確求解下列線性系統: Vjd+▽ψ(yj)=0, 或運用共軛梯度法找到滿足下述條件的方向dj:Vjdj+▽ψ(yj)≤min η, ▽ψ(yj)1+τ S2令αj=δmj, 其中mj是滿足下列式子的第一個非負整數m:ψ(yj+δmdj)≤ψ(yj)+μδm<▽ψ(yj), dj>S3 令yj+1=yj+αjdj

2 回歸模型及其優化算法

本節介紹三種統計學上重要的用于參數估計和變量預測的回歸模型，并且簡要回顧了近年來優化學者們對于求解相應模型的算法的研究進展.

2.1 線性回歸模型及其優化算法

在過去幾十年中，線性回歸模型作為統計學上的一個重要支柱，被廣泛應用于數據分析領域，特別是利用該模型預測連續的定量因變量. 但線性模型對收集到的樣本數據做了很多假設，并且該模型有可能產生不準確的預測變量.

2.1.1線性回歸模型和最小二乘估計

給定樣本(xi,yi)∈Rp×R,i=1, …,n,xi為p維自變量(特征向量或預測向量)，因變量(響應變量，觀測值)yi取常數值.線性回歸模型利用下述方程對觀測向量y=(y1,y2, …,yn)T∈Rn進行預測：

y=Xβ+

(4)

其中:X=(x1,x2, …,xn)T∈Rn×p是設計矩陣，β∈Rp是回歸系數，=(1,2, …,n)T∈Rn是隨機誤差.

2.1.2優化算法

因Lasso估計量有著優異的變量選擇和估計特性，并且Lasso問題即l1-正則最小二乘回歸問題是凸問題，便于求解. 近年來，學者們提出了許多求解Lasso估計量的優化算法，下面對歷年來提出的算法分為一階算法和二階算法分別進行概述.

Ⅰ) 一階算法. 2007年，Figueiredo等[89]提出了梯度投影法(GPSR)來求解l1-正則最小二乘回歸問題的等價形式—帶有界約束的二次規劃問題； 2008年，由于Lasso問題中的l1范數正則項具有可分結構，針對這一特殊結構, Wu等[90]提出貪婪坐標下降法； 之后于2010年，Friedman等[66]提出了循環坐標下降法. 2009年，Berg等[91]提出了譜投影梯度法(SPGL1)，SPGL1與GPSR方法都屬于梯度類算法： 通過將目標函數的梯度投影至某一可行集從而獲得相應迭代解； 同年，Wright等[92]提出了稀疏重構方法(SpaRSA). Wen等[93]提出了基于不動點連續方法(FPC)的改進方法-FPC_AS； Beck等[44]提出了快速軟閾值迭代算法(FISTA). 其中，FPC方法、 SpaRSA方法與FISTA方法的求解迭代形式相近，子問題都是求解目標函數于當前迭代步的二次近似展開式，詳細可見前面關于FISTA方法的具體描述. 2011年，Becker等[94]提出了NESTA算法，NESTA算法與FISTA都是基于Nestorov方法，但二者求解的目標問題形式不同并且算法生成的迭代序列也不同. Becker將NESTA算法與SPGL1、 SpaRSA、 FISTA、 FPC_AS等算法進行比較，實驗結果證明了SPGL1和NESTA算法具有更優異的數值表現，并且當需要求解高精度解時，NESTA算法更魯棒.

Ⅱ) 二階算法. 2001年，Chen等[95]將l1-正則最小二乘回歸問題轉化為線性規劃并運用經典的求解線性規劃的方法—內點法(IPM)求解. 內點法是一種原始-對偶問題優化方法即對原問題和對偶問題都進行優化迭代的算法，該算法根據問題的KKT條件，迭代尋求滿足KKT條件的解. 2007年，相較于文獻[95]提出的內點法，Kim等[96]設計了求解大規模l1-正則最小二乘問題的內點法，該方法在求解牛頓線性系統時采用的是預條件共軛梯度法，而在文獻[95]中采用共軛梯度法，并且兩者在迭代方法中采用的截斷準則也不一樣.數值實驗結果表明文獻[96]提出的算法更有優勢. 2014年，Fountoulakis等[97]對內點法進行了改進，求解過程中避免了對矩陣的直接計算和儲存，提出了無矩陣內點法mfIPM(matrix-free IPM). 2014年，Milzarek[98]提出了帶有多維度過濾全局化的半光滑牛頓法SNF. 2016年，Byrd等[99]提出了塊積極集算法BAS，在特定條件下該算法可看作與半光滑牛頓法相同. 2016年，Xiao等[100]提出了前向后向分裂-牛頓算法(FBS-Newton)求解Lasso問題，并證明了在一定條件下，FBS-Newton算法的收斂率可達二次. 2018年，Li等[80]提出了半光滑牛頓增廣拉格朗日方法求解Lasso問題，通過實驗證明了該方法遠遠優于FPC_AS、 ADMM、 APG、 mfIPM等先進算法.

2.2 廣義線性模型及其優化算法

2.2.1廣義線性模型

ⅰ) 廣義線性模型. 前面介紹的線性回歸模型適用于因變量是連續的定量變量的情況，然而在實際應用中存在著因變量是二值變量或計數變量即因變量yi取值是0或1，或者只取整數1，…，m中的一個值，此時線性回歸模型不再適用于刻畫變量xi和yi之間的關系，許多學者對線性回歸模型進行推廣得到廣義線性模型. 廣義線性模型的概念首先由Nelder等[101]提出，之后Mccullagh等[102]將其發展成系統的理論. 廣義線性模型在生物、 醫學和經濟等領域都有重要的應用意義，其具體形式如下：

2)y1, …,yn相互獨立，且yi的分布f(yi,μi,φ)是指數型分布.

關于指數型分布的具體定義，可參考文獻[7, 103]. 常見的重要分布如正態分布、 二項分布、 泊松分布均屬于指數型分布.

常見的聯系函數有如下幾種：

② Probit函數：g(p)=Φ-1(p)，Φ為N(0, 1)的分布函數;

③ ln-ln函數:g(p)=ln(-ln(1-p)).

2.2.2優化算法

近年來，隨著學者對序列二次近似方法(SQA)/鄰近牛頓法的改進，使得求解l1-正則邏輯回歸模型問題的方法不斷得到改善，如Yue等[73]提出了IRPN算法并比較了IRPN算法與FISTA算法[44]、 SpaRSA算法[92]、 坐標梯度下降算法CGD[119]、 newGLMNET[67]的數值表現，結果表明IRPN算法是更高效的，其可以在更短的時間內得到目標精度的解. 之后，Mordukhovich等[74]在理論和實驗上進一步對鄰近牛頓法改進，最后通過數值實驗證明了改進后的鄰近牛頓法比IRPN算法更高效求解l1-正則邏輯回歸模型問題.

2.3 分位數回歸模型及其優化算法

2.3.1分位數回歸模型

分位數回歸模型首先由Koenker等[120]于1978年提出. 分位數回歸模型是統計學和經濟計量學上一種應用廣泛的數據分析工具. 分位數回歸模型應用于許多研究領域如經濟[121- 122]、 生存分析[123- 124]、 微陣列研究[125]、 成長曲線圖[126-127]等.

分位數回歸模型不僅能夠度量解釋變量(自變量，回歸變量)對響應變量(因變量)的影響，并且能夠度量解釋變量對因變量分布上尾和下尾的影響. 其主要有如下兩個優點： Ⅰ)相較于線性回歸模型中的最小二乘估計，分位數回歸估計對線性回歸模型中的隨機誤差i的分布不需作任何假定要求，當隨機誤差i不服從高斯分布時，分位數回歸估計仍具有更優良的統計特性[120]； Ⅱ)分位數回歸是對所有分位數進行回歸，因此對數據中出現的異常點具有魯棒性.

① 分位數回歸模型. 設Y為實值隨機變量，分布函數為F(y)=P(Y≤y)，對任意的0<τ<1，Y的τ分位數為QY(τ)=inf {y:F(y)≥τ}.Koenker和Basset得到下述結果：

其中:ρτ(u):=(τ-I(u≤0))u是損失函數；I是指示函數.給定X和條件累積分布函數F(y|X)，則條件τ分位數QY(τ|x)=inf {y:F(y|X=x)≥τ}，并有：

考慮線性分位數回歸即QY(τ|x)是x的線性函數，不失一般性，假設QY(τ|x)=xTβ*，可得下列關于β*的分位數估計：

當τ=0.5時，F-1(τ)是中位數，此時的分位數回歸稱做中位數回歸或最小一乘回歸(最小絕對離差回歸LAD).

② 分位數回歸模型的研究現狀. Li等[128]，Wu等[129]首先研究了l1-正則分位數回歸模型，并且Wu等[129]證明了在較弱的條件下，帶有非凸正則項SCAD的分位數回歸模型具有漸近Oracle性質. 2012年，Wang等[130]證明了在一定條件下，即使在隨機誤差i服從重尾分布時，帶有非凸正則項(SCAD和MCP)的分位數回歸模型仍具有Oracle性質. 2013年，Wang[131]研究證明了l1-正則最小絕對離差回歸估計量在大概率下具有近似Oracle性質. 2014年，Fan等[132]研究了加權l1-正則分位數回歸模型，并將該模型用于處理重尾分布的高維數據以消除偏差，證明了在對隨機誤差i極弱的假設條件下，該模型具有Oracle性質及相應估計量具有漸近正態性.

2008年，Zou等[133]提出了復合分位數回歸模型并進一步研究了帶有自適應Lasso正則項的復合分位數回歸模型. 2010年，Kai等[134]將復合分位數回歸應用到非參數回歸中，相較于局部多項式回歸，提出了一種新的非參數回歸工具： 局部復合分位數回歸光滑化方法，并證明了針對誤差滿足不同非正態分布的情況下，該方法比局部多項式回歸方法更有效. 2011年，Kai等[135]進一步將復合分位數回歸應用于半參數回歸模型中. 2016年，Li等[136]研究了在Harris回歸馬爾可夫鏈框架下，針對非線性和非參數模型的局部多項式復合分位數回歸方法.

2.3.2優化算法

相較于最小二乘回歸和邏輯回歸問題，由于分位數回歸模型中的損失函數是非光滑且不可分的，求解該模型是更具有挑戰性的. 一種標準的求解方法是將分位數回歸問題轉化為線性規劃問題，然后運用優化軟件包求解變形后的線性規劃問題. 如2005年，Koenker等[137]應用內點法求解l1-正則分位數回歸問題，但內點法只適用于小規模或中規模問題. 受到Efron提出的求解Lasso問題的解路徑的最小角回歸算法[138](LARS算法不同于一般的優化迭代算法，該算法從空集開始逐步構造最優稀疏解，可看作一種序貫方法)的啟發，Li等[128]提出一種類似LARS的算法來有效求解l1-正則分位數回歸問題的解路徑: 分片線性路徑. 基于Zou等[139]提出的局部線性近似方法(LLA)，Fan等[140]提出了求解正則項是折疊凹性懲罰函數的正則化分位數回歸模型的方法: 兩步局部線性近似方法，并證明在大概率下，Oracle估計量可以在LLA算法中的兩步內獲得. 2015年，Peng等[141]提出了求解帶有非凸正則項(SCAD和MCP)的分位數回歸模型問題的快速迭代坐標下降方法(QICD)并證明了該算法產生的任何子序列收斂于一個穩定點，但理論上并不能保證該算法產生的解是正確的解. 2017年，Yi等[142]對帶有彈性網絡正則項的分位數回歸模型提出了半光滑牛頓坐標下降方法. 2018年，Gu等[58]提出了兩種基于交替方向乘子法的算法來求解帶有權重彈性網絡正則項的分位數回歸模型問題： 稀疏坐標下降ADMM算法(scdADMM)和半鄰近交替方向乘子法(SPADMM)，并發現當特征維數p大時scdADMM算法比SPADMM算法更高效，而p較小時SPADMM算法更有效. 2019年，Zhang等[143]提出了多階段凸松弛鄰近對偶半光滑牛頓法(MSCRA_PPA： 凸松弛和鄰近對偶半光滑牛頓法相結合，凸松弛是現今求解非凸問題的一種常見方法)求解高維下正則項是零范數的正則化分位數回歸估計量，并將MSCRA_PPA與多階段凸松弛內點法(MSCRA_IPM)、 多階段凸松弛交替方向乘子法(MSCRA_ADMM)進行比較，數值實驗證明了MSCRA_PPA可以花近乎其他算法一半的時間來得到同樣質量的解.

3 研究展望

在大數據時代，正則化回歸模型作為高維數據分析的重要工具可以用于變量選擇和參數估計，高效處理高維數據需要同時提高統計回歸模型的精確性和優化算法的計算效率. 本研究對統計學上的回歸模型和相關優化算法進行了簡要的概述. 主要介紹了三種重要且應用廣泛的回歸模型，并回顧了以l1范數為主的凸正則項和以折疊凹性懲罰函數(SCAD和MCP)為主的凹正則項.

以往，因傳統的連續優化理論與算法的局限性，優化學者更多集中于研究求解目標函數具有凸性和光滑性的回歸模型問題，如帶有凸正則項的最小二乘回歸問題和邏輯回歸問題. 接下來，對未來研究方向提出一些展望.

1) 相較于線性回歸模型，分位數回歸模型更適合于處理高維數據，因高維數據往往具有重尾分布等特征，而此時線性回歸模型對偏差分布需服從正態分布的要求無法得到滿足. 因此，對于求解分位數回歸模型問題的高效算法需要更多的研究.

2) 實際應用中，無論是對于回歸模型中的損失函數還是正則項，函數是非凸的情況是很普遍的，如機器學習中的經驗風險最小化中的損失函數往往是非凸函數. 近年來作為穩健處理高維數據的一種有效方法——低秩矩陣恢復，其中非凸函數可以更好近似秩函數. 深度學習算法中大多需要通過求解非凸優化問題來獲得最優的神經網絡參數. 另外，相較于凸正則項，以折疊凹性懲罰函數為代表的凹正則項具有更好的統計特性，如具有變量選擇一致性以及可以消除Lasso模型內蘊的估計偏差. 然而，現今非凸優化理論需要進一步完善，如非凸優化算法全局收斂性的分析，大部分時候非凸優化算法只能得到局部解，所得到的解離全局最優解有多近以及算法的魯棒性等都需要進一步的研究.