張　乾，何　岸，何洪津

(杭州電子科技大學(xué)理學(xué)院，浙江 杭州 310018)

?

一種解Dantzig-Selector模型的快速分解算法

張乾，何岸，何洪津

(杭州電子科技大學(xué)理學(xué)院，浙江 杭州 310018)

基于增廣拉格朗日法提出了一種快速分解算法求解Dantzig-Selector模型.與經(jīng)典的乘子交替方向法相比，新算法的每個(gè)子問題都具有更簡單易行的迭代格式.通過測試兩種不同類型的隨機(jī)數(shù)據(jù)，相應(yīng)的數(shù)值計(jì)算結(jié)果表明，算法在CPU運(yùn)行時(shí)間方面有較明顯的優(yōu)勢.

Dantzig-Selector模型；增廣拉格朗日方法；乘子交替方向法；分解算法

0　引　言

線性回歸是一類非常經(jīng)典的數(shù)學(xué)模型，它在信號(hào)處理、機(jī)器學(xué)習(xí)以及統(tǒng)計(jì)學(xué)習(xí)中有著極其廣泛的應(yīng)用.由于壓縮感知理論[1]的提出，尋找欠定線性回歸模型的稀疏解成為近年來最熱門的研究課題之一.然而，直接尋找滿足線性方程組的稀疏解是一個(gè)NP-難問題.為此，研究者提出了一系列凸松弛優(yōu)化模型，例如基追蹤模型[2]和LASSO模型[3].2007年，Candes和Tao針對(duì)超欠定線性方程組問題進(jìn)一步提出了更穩(wěn)健的Dantzig-Selector模型[4].與LASSO模型相比，Dantzig-Selector模型結(jié)構(gòu)更為復(fù)雜，從而給模型的求解增加了很大的困難.如何設(shè)計(jì)簡單高效的算法求解該模型是極具挑戰(zhàn)性的研究課題之一.本文首先通過引入兩個(gè)新的輔助變量，將Dantzig-Selector模型等價(jià)轉(zhuǎn)化為弱分離的優(yōu)化問題.然后，基于增廣拉格朗日方法，充分利用新模型的結(jié)構(gòu)特點(diǎn)設(shè)計(jì)出一種快速、可行的分解算法，使得算法的子問題都具有顯式表達(dá)式.相比已有的分裂算法，本文算法在CPU運(yùn)行時(shí)間上有較明顯的優(yōu)勢.

1　問題描述

經(jīng)典的線性回歸問題可刻畫為：

y=Xβ+ε.

(1)

式中，β∈Rp是未知的回歸系數(shù)，X∈Rn×p是一個(gè)有界的線性算子(或設(shè)計(jì)矩陣)，y∈Rn是帶噪音的觀測向量，ε是服從正態(tài)分布N(0,σ2I)的可加噪音.為了克服式(1)的病態(tài)性給模型求解過程帶來的不穩(wěn)定性，通常在最小二乘模型上引入Tikhonov正則化技術(shù).但是，當(dāng)β具有稀疏結(jié)構(gòu)時(shí)，且當(dāng)式(1)是一個(gè)超欠定線性系統(tǒng)，即n?p時(shí)，傳統(tǒng)的最小二乘模型及Tikhonov正則化得到的解無法滿足需求.2007年，Candes和Tao在LASSO模型的基礎(chǔ)上提出了Dantzig-Selector模型[4]：

(3)

上述模型是標(biāo)準(zhǔn)的帶線性約束的凸優(yōu)化模型，而增廣拉格朗日函數(shù)是求解此類問題的經(jīng)典方法之一.但是，直接使用增廣拉格朗日方法會(huì)使得子問題變量耦合在一起，不利于算法的實(shí)現(xiàn).因此，根據(jù)目標(biāo)函數(shù)的可分離性，文獻(xiàn)[6]利用乘子交替方向法進(jìn)行求解，且取得了良好的數(shù)值效果.美中不足之處是，直接利用乘子交替方向法求解，其關(guān)于β部分的子問題沒有顯式表達(dá)式，這給算法的實(shí)現(xiàn)增加了一定難度.隨后，文獻(xiàn)[7]巧妙地將與β有關(guān)的子問題進(jìn)行線性化，提高了乘子交替方向法的可操作性，也在一定程度上減少了算法的計(jì)算時(shí)間.但他們的方法要求估計(jì)一個(gè)矩陣的譜半徑，這也并非一件易事.基于上述原因，本文尋求一種新的簡單易行的算法.

2　變量分離算法

針對(duì)模型(3)，引入增廣拉格朗日函數(shù)：

(4)

式中，γ>0是罰參數(shù)，λ表示拉格朗日乘子.給定(βk,λk),文獻(xiàn)[5]中的乘子交替方向法的迭代格式可描述為：

(5)

(6)

λk+1=λk-γ(XTXβk+1-XTy-zk+1)

(7)

顯然，乘子交替方向法主要的任務(wù)就是求解子問題(5)和問題(6).由于子問題(5)可以顯式求解，子問題(6)的求解直接決定了算法的有效性.但由于XTX在通常情況下并不是一個(gè)單位矩陣，因此，式(6)并沒有顯式解，此時(shí)必須借助一些迭代算法進(jìn)行近似求解，這為算法的實(shí)現(xiàn)增加了較大的難度.基于乘子交替方向法的思想，下面提出一種新的分解算法，使得新算法所有子問題具有顯式迭代式.

為了簡化數(shù)學(xué)符號(hào)，記A∶=XTX和b∶=XTy.通過引入輔助變量x和z，模型(2)可轉(zhuǎn)化為弱分離優(yōu)化問題：

(8)

式中，ρ>0是一個(gè)懲罰因子.針對(duì)模型(8)，相應(yīng)的增廣拉格朗日函數(shù)為：

(9)

式中，λ和γ分別代表拉格朗日乘子和罰參數(shù).根據(jù)乘子交替方向法的Gauss-Seidel迭代思想，給定(xk,βk,λk)，分別對(duì)式(9)中的變量進(jìn)行交替求解，即得如下迭代格式：

(10)

(11)

(12)

λk+1=λk-γ(Aβk+1-b-zk+1)

(13)

類似迭代格式(5)～(7),上述方法的主要工作量集中在子問題(10)～(12).下文中，本文主要討論式(10)～(12)的具體表達(dá)式.

首先，關(guān)于變量z的子問題(10)歸結(jié)為：

(14)

(15)

式中，diag(D)表示由D矩陣對(duì)角元素構(gòu)成的列向量.

(16)

最后，將(zk+1,xk+1)代入(12)式，簡化關(guān)于β的子問題為：

(17)

(18)

由上述分析可見，本文新提出的算法其每個(gè)子問題都能顯式求解，這相比文獻(xiàn)[5]中的乘子交替方向法更容易實(shí)現(xiàn).當(dāng)矩陣A具有某些特殊結(jié)構(gòu)時(shí)，可以借助一些快速矩陣求逆方法對(duì)式(18)進(jìn)行求解(例如：快速傅立葉變換).當(dāng)矩陣A規(guī)模較大且無特殊結(jié)構(gòu)時(shí)，(γATA+ρI)的逆矩陣求解則會(huì)比較耗時(shí).此時(shí)，建議采用快速、成熟的共軛梯度法[8]進(jìn)行近似求解，即：

βk+1=βk-αkgk

(19)

當(dāng)處理(18)的病態(tài)情形時(shí)，在共軛梯度法中引入預(yù)處理迭代技術(shù)增加算法的穩(wěn)定性和可靠性，這也是新算法的一個(gè)潛在優(yōu)勢.綜上分析，本文的分解算法過程可描述如下：

新的快速分解算法解Dantzig-Selector模型的算法步驟如下：

1) 輸入初始迭代點(diǎn)(x0,β0,λ0)，選擇參數(shù)δ,ρ,γ>0；

2) 通過式(15)更新zk+1；

3) 通過式(16)更新xk+1；

4) 通過式(18)或(19)更新βk+1；

5) 通過式(13)更新λk+1.

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

本節(jié)對(duì)新提出的算法進(jìn)行數(shù)值模擬來檢驗(yàn)算法的優(yōu)劣.考慮到β子問題的計(jì)算方式，新算法分為式(18)的精確計(jì)算(簡記為“New-accurate”)和式(19)的一步迭代近似計(jì)算(簡記為“New-one”).同時(shí)，將新算法與文獻(xiàn)[5]中的乘子交替方向法(簡記為“ADMM”)進(jìn)行比較.所有算法都用MATLAB進(jìn)行編程實(shí)現(xiàn)，其中文獻(xiàn)[5]的ADMM算法程序由該文作者提供.

類似文獻(xiàn)[4]的測試情況，本節(jié)仍考慮兩種類型的系數(shù)矩陣X：單位列向量系數(shù)矩陣和行正交系數(shù)矩陣.實(shí)驗(yàn)中設(shè)定初始值迭代點(diǎn)為(x0,β0,λ0)=(0,0,0).新算法(New-one和New-accurate)迭代的停止條件設(shè)定為：

表1和表2的數(shù)值結(jié)果表明New-one算法和ADMM在得到質(zhì)量相當(dāng)?shù)慕平鈺r(shí)，New-one所需的CPU運(yùn)行時(shí)間比后者要少，這進(jìn)一步從數(shù)值上驗(yàn)證了本文新方法在求解子問題方面的優(yōu)勢.若直接使用式(18)更新βk+1，則需要計(jì)算矩陣的逆.一般情況下，計(jì)算矩陣的逆代價(jià)比較高，從New-accurate算法的數(shù)值結(jié)果也表明其CPU運(yùn)行時(shí)間較長.同時(shí)，本文處理的問題非常病態(tài)，故矩陣直接求逆的迭代方式導(dǎo)致計(jì)算結(jié)果不如另外兩種方式，因而進(jìn)一步說明式(19)的重要性.但是，在矩陣維數(shù)比較低時(shí)，New-accurate的近似解效果還是令人滿意，從某種程度上說，當(dāng)矩陣X具有某種特殊結(jié)構(gòu)時(shí)，或許可采用快速傅立葉變換實(shí)現(xiàn)高精度的快速求解.

表1　當(dāng)X為單位列向量矩陣時(shí)，3種算法數(shù)值結(jié)果

表2　當(dāng)X為行正交矩陣時(shí)，3種算法數(shù)值結(jié)果

圖1　X為行正交矩陣時(shí)的迭代次數(shù)比較

圖2　X為單位列向量矩陣時(shí)的迭代次數(shù)比較

圖3　X為行正交矩陣時(shí)的計(jì)算誤差比較

圖4　X為單位列向量矩陣時(shí)的計(jì)算誤差的比較

圖5和圖6分別畫出了當(dāng)(n,p,s)=(72,256,8)時(shí)兩種不同類型系數(shù)矩陣對(duì)應(yīng)的數(shù)值解情況.兩幅圖可看出New-one恢復(fù)出來的數(shù)值解令人滿意.

圖5　X為行正交矩陣時(shí)的近似解

圖6　X為單位列向量矩陣時(shí)的近似解

4　結(jié)束語

鑒于直接使用乘子交替方向法求解Dantzig-Selector模型時(shí)子問題不具有顯式迭代格式的缺陷，本文引入一種新的模型刻畫技巧，提出了一種子問題簡單易行的分解算法.通過測試兩種不同類型的隨機(jī)數(shù)據(jù)，初步的數(shù)值結(jié)果表明新算法在CPU運(yùn)行時(shí)間方面有較明顯的優(yōu)勢.

[1]DONOHO D L.Compressed sensing[J]. Information Theory, IEEE Transactions on,2006,52(4):1289-1306.

[2]CHEN S S, DONOHO D L, SAUNDERS M A. Automatic decomposition by basis pursuit[J]. SIAM Journal on Scientific Computing,1998,20(1):33-61.

[3]TIBSHIRANI R. Regression shrinkage and selection via the lasso[J]. Journal of the Royal Statistical Society B,1996,58(1):267-288.

[4]CANDES E, TAO T. The Dantzig selector: statistical estimation when p is much larger than n[J]. Annals of Statistics,2007,35(6):2313-2351.

[5]LU Z S, PONG T K, ZHANG Y. An alternating direction method for finding Dantzig selectors[J]. Computational Statistics & Data Analysis,2012,56(12):4037-4046.

[6]GABAY D, MERCIER B. A dual algorithm for the solution of nonlinear variational problems via finite element approximations[J]. Computers & Mathematics with Applications,1976,2(1):16-40.

[7]WANG X, YUAN X. The linearized alternating direction method of multipliers for Dantzig selector[J]. SIAM Journal on Scientific Computing,2012,34(5):A2792-A2811.

[8]袁亞湘,孫文瑜.最優(yōu)化理論與方法[M].北京:科學(xué)出版社,1996:186-194.

A Fast Decomposition Algorithm for Finding Dantzig-Selectors

ZHANG Qian, HE An, HE Hongjin

(SchoolofScience,HangzhouDianziUniversity,HangzhouZhejiang310018,China)

Based on the augmented Lagrangian method, this paper introduces a fast decomposition algorithm for solving the Dantzig selector model. Comparing with the classical alternating direction method of multipliers, all subproblems of the proposed algorithm have closed-form solutions so that the new algorithm is easily implementable in practice. Finally some preliminary numerical results show that our new algorithm has superiority in terms of taking less CPU time by testing two types of synthetic data sets.

Dantzig-Selector model; augmented Lagrangian method; alternating direction method of multipliers; decomposition method

10.13954/j.cnki.hdu.2016.01.020

2015-06-22

浙江省自然科學(xué)基金重點(diǎn)資助項(xiàng)目(LZ14A010003)；浙江省大學(xué)生新苗人才計(jì)劃資助項(xiàng)目(2015R407038)

張乾(1993-)，女，河南周口人，在讀本科生，數(shù)值優(yōu)化和變分不等式.通信作者：何洪津講師，E-mail: hehjmath@hdu.edu.cn.

O221.2

A

1001-9146(2016)01-0097-06