1、337Journalof BeijNormalUniversity(NaturalScience)59(2)2023-04北京师范大自然科学版）高维惩罚分位数回归及优化算法袁攀旭1)罗敬宣1)岳莉莉2)李高荣1)+（1)北京师范大学统计学院，10 0 8 7 5，北京；2)南京审计大学统计与数据科学学院，2 118 15，江苏南京）摘要针对具有异常值或离群点的高维数据线性回归模型，提出了一种基于误差函数正则化的惩罚分位数回归的新方法，与经典的L惩罚方法相比，新方法具有更好的稳健性以及更小的估计偏差和预测误差；为解决分位数损失函数非光滑性与误差函数非凸性所带来的计算挑战，结合迭代再加权L算法以及
2、ADMM算法，提出了一种有效的IRW-ADMM算法，并对回归系数进行了求解模拟结果表明，与已有的惩罚分位数回归方法相比，新方法在参数估计和变量选择等方面均具有更好的表现.将新方法应用于核黄素基因数据分析，以证实其有效性和可行性。关键词高高维数据；分位数回归；变量选择；稳健估计；误差函数；IRW-ADMM算法中图分类号0212.1D01:10.12202/j.0476-0301.20223110引言现代计算机技术的飞速发展，极大提升了各行各业的数据收集能力，产生了海量的大数据；数据的特征不仅体现为数据量的增大，更表现为数据维度的增高，从而导致收集到的数据中往往包含众多异常值或离群点.为更好地探索
3、在异常值或离群点存在的情况下响应变量与解释变量之间的相关关系，Koenker等提出了分位数回归方法，被证实是一种有效的稳健估计方法，并被广泛应用于基因组学、经济学、金融学和社会科学等领域，截至目前已取得丰富的研究成果 2-4,然而，在当今大数据时代，传统的分位数回归无法直接用来处理高维数据，这对统计学者提出了新的要求和挑战，呕需发展适用于高维数据的稳健估计方法。庆幸的是，尽管收集到的数据维度很高，但在多数情况下仅有少部分协变量与响应变量之间存在显著性统计关系，即稀疏性假设.基于该假设，在过去20多年间，统计学者们已经发展了一类惩罚最小二乘估计方法.该类方法可同时进行变量选择和参数估计，其在实际
4、应用中具有较满意的预测性能，已在众多科学领域得到了广泛应用，其中具有代表性的方法包括Tibshirani5提出的LASSO方法，Fan等 6 提出的SCAD方法，Zou等 7 提出的弹性网方法，Zoul8提出的自适应LASSO方法，Zhang9提出的MCP方法等.Fan等 10-1对这类方法做了系统的总结.在高维背景下，数据维度的增加不仅导致协变量之间存在伪相关，也使得惩罚最小二乘估计对数据进行拟合时，更容易受到数据中异常值或离群点的影响，从而导致参数估计和变量选择结果失真.因此，需要发展新的高维稳健估计方法处理复杂高维数据分析中异常值或离群点造成的影响。近年来，在高维回归模型下利用惩罚方法研
5、究分位数回归模型的参数估计和变量选择问题，已成为统计学领域的研究热点.例如：在固定维情形下，Wu等 12 研究了SCAD惩罚和自适应LASSO惩罚分位数估计的Oracle性质及变量选择相合性；Belloni等 13在高维情形下提出了L惩罚分位数回归方法，并证明了估计量的相合性；Wang等 14、Wangl15、Fa n 等 16 分别在不同情形下发展了惩罚分位数回归方法，并证明了估计量的Oracle性质和变量选择相合性.与最小二乘相比，分位数回归的损失函数是非光滑函数，这在高维情形下为惩罚分位数回归问题的求解提出了新的挑战.已有研究提出了几种优化算法求解惩罚分位数回归问题.例如：Koenker
6、等 17 基于内点法构造了适用于分位数回归和惩罚分位数回归模型的求解算法；Li等 18 利用线性规划对L惩罚分位数回归问题进行求解，并且给出了解的路径；Wu等 19基于坐标下降算法提出一种求解L惩罚分位数回归模型的贪婪算*国家自然科学基金资助项目（12 2 7 10 46，1197 10 0 1，12 0 0 12 7 7）；国家自然科学基金重点资助项目（12 1310 0 6）十通信作者：李高荣（197 6 一），男，教授，博导。研究方向：高维统计、非参数统计、复杂数据分析，以及统计学习等.E-mail：收稿日期：2 0 2 2-0 4-30338第59 卷北京师范大自然科学版）法；Peng
7、等 2 0 基于坐标下降算法给出了一般惩罚分位数回归模型的求解算法；Gu等 2 1 对求解高维惩罚分位数回归模型提出了近端ADMM算法和稀疏坐标下降ADMM算法改进的优化算法；Wang等 2 2 进一步将协变量之间的结构信息纳人分位数回归模型中，提出自适应分位数回归的ADMM求解算法.高维分位数回归模型的估计和计算等方面的研究目前仍是热点.为进一步扩展高维惩罚分位数回归的方法论，本文提出一种基于误差函数正则化的高维惩罚分位数回归方法.与L惩罚分位数回归方法相比，新提出的方法不仅导致更小的估计偏差，而且能够准确识别显著变量.此外，误差函数能够很好地逼近单位跳跃函数，故在一定程度上可以很好地近似L
8、.惩罚函数.然而，误差函数非凸性和分位数损失函数非光滑性在数值优化方面为所提方法提出了重大挑战.为解决这一困难，本研究通过结合Candes等 2 3的迭代再加权（iterativereweighted，IR W)L算法和ADMM算法，提出迭代再加权ADMM(IRW-ADMM)算法，以获得基于误差函数惩罚分位数估计的数值解.IRW-ADMM算法的核心步骤有3个：1)给定初值，计算权重，将优化问题转化为求解加权L惩罚分位数回归问题；2)利用ADMM算法得到加权L惩罚分位数回归的数值解；3)更新权重，迭代求解，直至收敛。模拟结果表明，该算法与已有惩罚分位数回归方法相比，在变量选择和参数估计等方面均具
9、有更好的表现，因此将IRW-ADMM算法用于核黄素基因数据分析，以识别与核黄素合成有关的基因.1基于误差函数的惩罚分位数回归分位数回归经常被用于研究一组协变量对响应变量条件分布的影响，不仅对异常值或离群点具有稳健性，而且适用于异方差的情形.对于响应变量YER和p维协变量X=(Xi,X,)TeR,令Fy=P(YyX=x)表示给定X=x的条件下y的条件分布，Qr(T|x)=inf(y:F(y|x)T)表示t分位数,其中T E(0,1).线性分位数回归模型假定Qr(ylx)=xT(T),其中(t)表示p维回归系数向量.令(xi,y),i=1,n)表示从总体(X,Y)中收集到的独立同分布观测样本。对高
10、维线性分位数回归模型，通常采用惩罚分位数回归方法得到回归系数(T)的稀疏估计，即求解优化问题minn=1j=1式中：p.(u)=u(t-I(u0为惩罚函数；入 0 为调节参数.常用的惩罚函数有如下3种。1)L,或 LASSO:p,(lul)=|ul.2)SCAD:对常数a2,定义为八ul,u入,2a|u|-ul-22p;(u)=入an.23）M C P：对常数a0,定义为ullula入,pa(u)=2a(a入)/2,ul a入.L惩罚函数是凸函数，可利用凸优化算法进行求解，但所得估计量有偏差.SCAD和MCP惩罚函数均为凹函数，可纠正L,惩罚函数对较大回归系数所造成的估计偏差，并在一定条件下可
11、证明所得估计量具有Oracle性质，但不具有光滑性.不同于已有惩罚函数，本文提出了一种新的基于误差函数的非凸正则化方法以获得惩罚分位数估计。考虑优化问题min(1)n=1式中:J。():=。(I,I)为误差函数正则项；。(x)=j=1e-r/dt表示误差函数；0为参数.本文将求解优化问题式(1)所获得的惩罚分位数估计，称为误差函数正则化(error function regularization,EFR)估计.函数。()与J。)具有如下4条性质.性质1误差函数。(x)处处可导,并且有。(x)=e-t2/o2性质2对任意的xER，误差函数。(x)满足c V1-e-ax。(x)c V1-e-b p
12、,其中a=1/,b=元/(42),C=QV元/2.性质3函数J.在R上为凹函数，即对任意的tE0,1和,e R,有J。(t+(1-t)tJ。()+(1-t)J.().性质4函数J。满足次可加性，即对任意的和ERP,有J。(+)J.()+J。(B).Guo等 2 4给出了性质1 4的详细证明.性质1说明误差函数具有光滑性，并且当x+8时其导函数(x)趋于0.性质2 表明误差函数在x+时有界.性质12 使得误差函数与L惩罚函数存在本质339袁攀旭等：高维惩罚分位数回归及优化算法第2 期区别.性质3表明误差函数在(0,+8)上为凹函数.性质4说明误差函数满足三角不等式.此外，误差函数还具有非常好的渐
13、近性质，当分别趋于0 或+时，误差函数正则化分别等价于L或L正则化。命题1对任意的非零向量ER，有如下2 个结论：1)当+8 时,有J();2)当0*时,有J。()/2-1元l.命题1的证明参见文献 2 4.命题1表明误差函数介于L惩罚函数和L惩罚函数之间，并且对于较小的，误差函数可以很好地近似L惩罚函数；而当+8时，误差函数等价于L惩罚函数.图1展示了当入=1时几种不同的惩罚函数.6LEFR:G=0.25EFR:=1.0EFR:G=2.04SCAD3P210-6-4-20246u图1不同的惩罚函数图1中，横坐标为自变量u,其取值范围为-6,6,纵坐标为惩罚函数p(u)的取值.为方便进行对比，
14、所有函数的图形进行了适当的缩放，使其均经过（1，1)点.由图1可见，对于较小的，误差函数能够很好地近似L惩罚函数，并且与其他惩罚函数相比，其“偏差”更小；随着的增大，误差函数逐渐向L惩罚函数“靠近”.2IRW-ADMM算法分位数损失函数非光滑性和误差函数非凸性，对优化问题式(1)的数值求解提出了极大挑战.为有效求解优化问题式（1），本文结合迭代再加权L算法和ADMM算法，提出IRW-ADMM算法对优化问题式（1)进行求解。由性质1可知，(x)=e-r/o.因此可通过迭代再加权惩罚L算法求解优化问题式（1），按如下3个步骤进行步骤1对T E(O，1),给定(T)的初始值和90.步骤2对11:1)
15、计算权重w=exp(-(9/o),j=1,p;2)求解如下加权L惩罚分位数回归问题min(2)=1j=1得到估计BO.步骤3若-(-1)1,则终止算法，输出,否则返回步骤2.式中：IIl表示向量的I,范数(q1);表示送代误差的阅值.由于缺少显式解，精确求解子问题式（2)是不切实际的。为解决这一困难，本文采用ADMM求解算法。即考虑优化问题min(3)=1j=1令Q(z)=n-Zp.(c),其中 z=-.y=(J1,yn)T,X=(xi,xn)T，因此，式（3)等价于带约束的优化问题min(Q.(z)+Allw oll),zs.t.:X+z=y,（4)式中aol=/Bgl,对固定的0.式(4)
16、的增广=1Lagrange 函数为SL(,z,0):=Q.(z)+w oll-(0,X+z-y)+号lX+z-yll2,式中一，ST-1TProx,.5,a=o,T一，T-1T-1一S故此，式（6）可表示为2(k+1)=Prox,.(y;-x+,no,i=l,.,n.与z步不同，迭代步骤（5)中的步无法推导出简单的显式解,因此，本文考虑在步的目标函数中添加近端（proximal)项，得到增广的步，其计算步骤为(k+1):=argmin(llw o ll-(0(),X)+(7)2式中：lvll：=0,m a x o u t e r O,m a x in n e r 0;1初始化外循环：1=1,给定初值(B(%)，2(0)，0);2当1i时,计算3w(=exp(-(0/);4初始化内循环：k=1,()o;5当k2时,更新6(+)(shrinks)+()+6y-X0)-2),),1jp7z(+1)(Prox,y.-(+1+-0l,nl),1in8(+1)()-y8(X(k+1)+z(k1)-y);9kk+l;10停止11(+)=(),I+1;12停止输出：0.3模拟研究为验证基于式(1)得到