高维统计学习笔记1——LASSO和Oracle性质

主要参考资料：Sara Van De Geer《Estimation and Testing Under Sparsity》

前言

当年Tibshirani提出LASSO的时候，未曾想到LASSO竟然成为了高维统计中一个非常重要的工具，这其中当然有一部分要归功于苏黎世联邦理工大学的美女教授Sara van de geer对LASSO理论的贡献。废话少说，开始学习。

1.高维统计的重要工具——极小化正则风险

当数据特征的数量ppp远大于我们所观测的样本量nnn时，我们称数据是高维的，如果把总体的特征看作参数β\betaβ，对参数进行估计和检验的一个有效方法是极小化正则风险。

Notation

损失函数 R:B→R,B⊂RpR:\mathcal{B}\rightarrow R, \mathcal{B}\subset\mathbb{R}^pR:B→R,B⊂Rp.
目标参数 β0=arg⁡min⁡β∈BR(β)\beta^0=\arg\min_{\beta\in\mathcal{B}}R(\beta)β0=argminβ∈BR(β).
经验损失函数 Rn:B→RR_n:\mathcal{B}\rightarrow RRn:B→R, based on n data points X1,X2,...,XnX_1,X_2,...,X_nX1,X2,...,Xn with n<pn<pn<p
惩罚 pen:B→[0,+∞)pen:\mathcal{B}\rightarrow[0,+\infty)pen:B→[0,+∞),实际上它是一个给定稀疏性后的惩罚，它对应着Rp\mathbb{R}^pRp上的一个范数。

如何估计参数β\betaβ呢？我们考虑的极小化正则风险：
β^:=arg⁡min⁡β∈B{Rn(β)+pen(β)}\hat{\beta}:=\arg\min_{\beta\in\mathcal{B}}\{R_n(\beta)+pen(\beta)\}β^:=argβ∈Bmin{Rn(β)+pen(β)}Van de geer 这本书的目的就正是研究这个估计，希望能有在很高的概率下有下面这个不等式：R(β^)≤R(β0)+Remainder，R(\hat{\beta})\leq R(\beta_0)+Remainder，R(β^)≤R(β0)+Remainder， RemainderRemainderRemainder是一个很小的数，它取决于β0\beta_0β0有多稀疏。

而她实际展示的结果是下面这个更一般的不等式：R(β^)≤R(β)+Remainder(β)，∀β.R(\hat{\beta})\leq R(\beta)+Remainder(\beta)，\forall\beta.R(β^)≤R(β)+Remainder(β)，∀β.这里Remainder(β)Remainder(\beta)Remainder(β)取决于β\betaβ的非零项的数目，其实这就是一个sharp oracle inequality，而一个non-sharp oracle inequality有下面的形式:R(β^)≤R(β0)+(1+η)(R(β)−R(β0))+Remainder(β),∀β.R(\hat{\beta})\leq R(\beta_0)+(1+\eta)(R(\beta)-R(\beta_0))+Remainder(\beta),\forall \beta.R(β^)≤R(β0)+(1+η)(R(β)−R(β0))+Remainder(β),∀β. 这里η>0\eta>0η>0是某个很小的正常数。

2.线性模型和LASSO

Notation

∣∣.∣∣n2=∣∣.∣∣22/n,||.||_n^2=||.||_2^2/n,∣∣.∣∣n2=∣∣.∣∣22/n,βS,j=βj1{j∈S},\beta_{S,j}=\beta_j1\{j\in S\},βS,j=βj1{j∈S},βS=(βS,1,...,βS,p)T,\beta_S=(\beta_{S,1},...,\beta_{S,p})^T,βS=(βS,1,...,βS,p)T,β−S=βSc.\beta_{-S}=\beta_{S^c}.β−S=βSc.
我们考虑模型Y=f0+ϵ.Y=f^0+\epsilon.Y=f0+ϵ.这里ϵ\epsilonϵ是一个期望为0的干扰，Y∈RnY\in\mathbb{R}^nY∈Rn。假设f0=Xβ0f^0=X\beta^0f0=Xβ0，那么这就是一个线性模型。假设ϵi\epsilon_iϵi之间不相关且方差为σ2\sigma^2σ2，Xn×pX_{n\times p}Xn×p是一个非随机的设计阵，这时我们注意到 E∣∣Ynew−Xβ^∣∣n2=E∣∣Xβ0−Xβ^∣∣n2+σ2.E|| Y_{new}-X\hat{\beta}||_n^2=E||X\beta_0-X\hat{\beta}||_n^2+\sigma^2.E∣∣Ynew−Xβ^∣∣n2=E∣∣Xβ0−Xβ^∣∣n2+σ2.
显然，我们希望∣∣Xβ0−Xβ^∣∣n2=op(1)|| X\beta_0-X\hat{\beta}||_n^2=o_p(1)∣∣Xβ0−Xβ^∣∣n2=op(1)。

如果n>pn>pn>p，最小二乘估计β^LS\hat{\beta}^{LS}β^LS是对参数β\betaβ的一个合理的估计，我们很容易得到E∣∣Xβ0−Xβ^LS∣∣n2=pnσ2,E||X\beta_0-X\hat{\beta}^{LS}||_n^2=\frac{p}{n}\sigma^2,E∣∣Xβ0−Xβ^LS∣∣n2=npσ2,当p≫np\gg np≫n且XXX行满秩时，有
E∣∣Xβ0−Xβ^LS∣∣n2=nnσ2=σ2≠o(1).E||X\beta_0-X\hat{\beta}^{LS}||_n^2=\frac{n}{n}\sigma^2=\sigma^2 =\not o(1).E∣∣Xβ0−Xβ^LS∣∣n2=nnσ2=σ2≠o(1).
因此在高维条件下最小二乘估计是不合适的。LASSO一个很好的性质是它能够产生稀疏的估计，这源于L1L_1L1约束的几何形状，这点想必大家都清楚，就不多说了。LASSO估计的定义是：β^=arg⁡min⁡β∈B{∣∣Y−Xβ∣∣n2+2λ∣∣β∣∣1}.\hat{\beta}=\arg\min_{\beta\in\mathcal{B}}\{||Y-X\beta||_n^2+2\lambda||\beta||_1\}.β^=argβ∈Bmin{∣∣Y−Xβ∣∣n2+2λ∣∣β∣∣1}.同样，我们去计算一下∣∣Xβ0−Xβ^∣∣n2||X\beta_0-X\hat{\beta}||_n^2∣∣Xβ0−Xβ^∣∣n2，我们希望它在高维的情况下依概率收敛到0，并且想知道速度有多快。

首先，由LASSO的定义显然有
∣∣Y−Xβ^∣∣n2+2λ∣∣β^∣∣1≤∣∣Y−Xβ∣∣n2+2λ∣∣β∣∣1,∀β||Y-X\hat{\beta}||_n^2+2\lambda|\hat{|\beta}||_1\leq||Y-X\beta||_n^2+2\lambda||\beta||_1,\forall \beta∣∣Y−Xβ^∣∣n2+2λ∣∣β^∣∣1≤∣∣Y−Xβ∣∣n2+2λ∣∣β∣∣1,∀β做一点简单的计算可得∣∣Xβ0−Xβ^∣∣n2≤∣∣Xβ0−Xβ∣∣n2+2ϵTX(β^−β)n+2λ(∣∣β∣∣1−∣∣β^∣∣1),||X\beta^0 -X\hat{\beta}||_n^2\leq||X\beta^0-X\beta||_n^2+\frac{2\epsilon^TX(\hat{\beta}-\beta)}{n}+2\lambda(||\beta||_1-||\hat{\beta}||_1),∣∣Xβ0−Xβ^∣∣n2≤∣∣Xβ0−Xβ∣∣n2+n2ϵTX(β^−β)+2λ(∣∣β∣∣1−∣∣β^∣∣1),为了去bound左边这个prediction error，我们需要bound后面三项，第一项暂且不用去管它，对第二项有
2ϵTX(β^−β)n≤2∣∣ϵTXn∣∣∞∣∣β^−β∣∣1,\frac{2\epsilon^TX(\hat{\beta}-\beta)}{n}\leq2||\frac{\epsilon^TX}{n}||_{\infty}||\hat{\beta}-\beta||_1,n2ϵTX(β^−β)≤2∣∣nϵTX∣∣∞∣∣β^−β∣∣1,我们不妨令集合
F={w:2∣∣ϵ(w)TXn∣∣∞<2λ0},\mathcal{F}=\{w:2||\frac{\epsilon(w)^TX}{n}||_{\infty}<2\lambda_0\},F={w:2∣∣nϵ(w)TX∣∣∞<2λ0},如果ϵ∼Nn(0,σ2I)\epsilon\sim\mathcal{N}_n(0,\sigma^2I)ϵ∼Nn(0,σ2I)，而且我们的数据阵经过了标准化,i.e.,∣∣X(j)∣∣n2=1||X^{(j)}||_n^2=1∣∣X(j)∣∣n2=1，则ϵTX(j)/n∼N(0,σ2n)\epsilon^TX^{(j)}/n\sim\mathcal{N}(0,\frac{\sigma^2}{n})ϵTX(j)/n∼N(0,nσ2)，这时P[Fc]≤∑j=1pP[∣ϵTX(j)n∣≥λ0]≤pe−nλ022σ2=eln⁡(p)−nλ022σ2.\mathbb{P}[\mathcal{F}^c]\leq\sum_{j=1}^{p}\mathbb{P}[|\frac{\epsilon^TX^{(j)}}{n}|\geq\lambda_0]\leq pe^{-\frac{n\lambda_0^2}{2\sigma^2}}=e^{\ln(p)-\frac{n\lambda_0^2}{2\sigma^2}}.P[Fc]≤j=1∑pP[∣nϵTX(j)∣≥λ0]≤pe−2σ2nλ02=eln(p)−2σ2nλ02.令λ0=σ2log⁡p+t2n\lambda_0=\sigma\sqrt{\frac{2\log p+t^2}{n}}λ0=σn2logp+t2，则有P[F]≥1−e−t2/2\mathbb{P}[\mathcal{F}]\geq1-e^{-t^2/2}P[F]≥1−e−t2/2，也就是说如果λ0≍log⁡(p)n\lambda_0\asymp\sqrt{\frac{\log(p)}{n}}λ0≍nlog(p)，那么P[F]≈1,\mathbb{P}[\mathcal{F}]\approx1,P[F]≈1,即∣∣ϵ(w)TXn∣∣∞=Op(log⁡(p)n).||\frac{\epsilon(w)^TX}{n}||_{\infty}=O_p(\sqrt{\frac{\log(p)}{n}}).∣∣nϵ(w)TX∣∣∞=Op(nlog(p)).
所以在集合F\mathcal{F}F上，
∣∣Xβ0−Xβ^∣∣n2≤∣∣Xβ0−Xβ∣∣n2+2λ0∣∣β^−β∣∣1+2λ(∣∣β∣∣1−∣∣β^∣∣1),||X\beta^0 -X\hat{\beta}||_n^2\leq||X\beta^0-X\beta||_n^2+2\lambda_0||\hat{\beta}-\beta||_1+2\lambda(||\beta||_1-||\hat{\beta}||_1),∣∣Xβ0−Xβ^∣∣n2≤∣∣Xβ0−Xβ∣∣n2+2λ0∣∣β^−β∣∣1+2λ(∣∣β∣∣1−∣∣β^∣∣1),(1)(1)(1)令β=β0\beta=\beta^0β=β0，则有
(2)0.5∣∣Xβ0−Xβ^∣∣n2+(λ−λ0)∣∣β^∣∣1≤(λ+λ0)∣∣β0∣∣1,0.5||X\beta^0 -X\hat{\beta}||_n^2+(\lambda-\lambda_0)||\hat{\beta}||_1\leq(\lambda+\lambda_0)||\beta^0||_1, \tag20.5∣∣Xβ0−Xβ^∣∣n2+(λ−λ0)∣∣β^∣∣1≤(λ+λ0)∣∣β0∣∣1,(2)如果取λ0≤0.5λ\lambda_0\leq0.5\lambdaλ0≤0.5λ，则，
(3)∣∣Xβ0−Xβ^∣∣n2+λ∣∣β^∣∣1≤3λ∣∣β0∣∣1.||X\beta^0 -X\hat{\beta}||_n^2+\lambda||\hat{\beta}||_1\leq3\lambda||\beta^0||_1. \tag3∣∣Xβ0−Xβ^∣∣n2+λ∣∣β^∣∣1≤3λ∣∣β0∣∣1.(3)所以我们关心的为问题变成了∣∣β0∣∣1||\beta^0||_1∣∣β0∣∣1到底有多大。这时，我们需要在∣∣β0∣∣1||\beta^0||_1∣∣β0∣∣1和β0TΣ^β0{\beta^0}^T\hat{\Sigma}\beta^0β0TΣ^β0之间建立起联系，记Σ^=XTX/n\hat{\Sigma}=X^TX/nΣ^=XTX/n。Sara van de geer(2007) 那篇文章里提出了Compatibility constant，它的定义是，对一个常数L≥1L\geq1L≥1和指标集SSS,

ϕ^2(L,S):=min⁡{∣S∣∣∣XβS−Xβ−S∣∣n2:∣∣βS∣∣1=1,∣∣β−S∣∣1≤L},\hat{\phi}^2(L,S):=\min\{|S|||X\beta_S-X\beta_{-S}||_n^2:||\beta_S||_1=1,||\beta_{-S}||_1\leq L\},ϕ^2(L,S):=min{∣S∣∣∣XβS−Xβ−S∣∣n2:∣∣βS∣∣1=1,∣∣β−S∣∣1≤L}, LLL一般被称作“拉伸因子”，直观上看，C-constant其实就是一个凸包到另一个拉伸后的凸包的距离。对任意β∗\beta^*β∗，取S∗={j:βj∗≠0},δ∗=β∗∣∣β∗∣∣1S^*=\{j:\beta^*_j=\not0\},\delta^*=\frac{\beta^*}{||\beta^*||_1}S∗={j:βj∗≠0},δ∗=∣∣β∗∣∣1β∗，那么显然有∣∣β∗∣∣12≤∣S∗∣∣∣Xβ∗∣∣n2ϕ^2(1,S∗).||\beta^*||_1^2\leq\frac{|S^*|||X\beta^*||_n^2}{\hat{\phi}^2(1,S^*)}.∣∣β∗∣∣12≤ϕ^2(1,S∗)∣S∗∣∣∣Xβ∗∣∣n2.当然，这个S∗S^*S∗通常是不知道的，如果考虑的是任意集合SSS，显然我们需要再添加一个条件∣∣β−S∗∣∣1≤3∣∣βS∗∣∣1||\beta_{-S}^*||_1\leq3||\beta_{S}^*||_1∣∣β−S∗∣∣1≤3∣∣βS∗∣∣1(这里的3我给的相当随意)，这样只需要令δ∗=β∗∣∣βS∗∣∣1\delta^*=\frac{\beta^*}{||\beta_S^*||_1}δ∗=∣∣βS∗∣∣1β∗，就有∣∣βS∗∣∣12≤∣S∣∣∣XβS∗−Xβ−S∗∣∣n2ϕ^2(3,S).||\beta_{S}^*||_1^2\leq\frac{|S|||X\beta_S^*-X\beta_{-S}^*||_n^2}{\hat{\phi}^2(3,S)}.∣∣βS∗∣∣12≤ϕ^2(3,S)∣S∣∣∣XβS∗−Xβ−S∗∣∣n2.有了C-Constant的概念，我们取S0={j:βj0≠0}S^0=\{j:\beta^0_j=\not0\}S0={j:βj0≠0}，对(1)做一点简单的变换，在λ>2λ0\lambda>2\lambda_0λ>2λ0的条件下可得，(4)∣∣Xβ0−Xβ^∣∣n2+λ∣∣β^−S0∣∣1≤3λ∣∣β^S0−βS00∣∣1,||X\beta^0 -X\hat{\beta}||_n^2+\lambda||\hat{\beta}_{-S^0}||_1\leq3\lambda||\hat{\beta}_{S^0}-\beta^0_{S^0}||_1,\tag4 ∣∣Xβ0−Xβ^∣∣n2+λ∣∣β^−S0∣∣1≤3λ∣∣β^S0−βS00∣∣1,(4)取δ=β^−β0∣∣β^S0−βS00∣∣1\delta=\frac{\hat{\beta}-\beta^0}{||\hat{\beta}_{S^0}-\beta^0_{S^0}||_1}δ=∣∣β^S0−βS00∣∣1β^−β0，可以得到∣∣β^S0−βS00∣∣12≤∣S0∣∣∣Xβ0−Xβ^∣∣n2ϕ^2(3,S0).||\hat{\beta}_{S^0}-\beta^0_{S^0}||_1^2\leq\frac{|S^0|||X\beta^0 -X\hat{\beta}||_n^2}{\hat{\phi}^2(3,S^0)}.∣∣β^S0−βS00∣∣12≤ϕ^2(3,S0)∣S0∣∣∣Xβ0−Xβ^∣∣n2.由基本不等式，
(5)4λ∣∣β^S0−βS00∣∣1≤8λ2∣S0∣ϕ^2(3,S0)+0.5∣∣Xβ0−Xβ^∣∣n2,4\lambda||\hat{\beta}_{S^0}-\beta^0_{S^0}||_1\leq\frac{8\lambda^2|S^0|}{\hat{\phi}^2(3,S^0)}+0.5||X\beta^0 -X\hat{\beta}||_n^2,\tag54λ∣∣β^S0−βS00∣∣1≤ϕ^2(3,S0)8λ2∣S0∣+0.5∣∣Xβ0−Xβ^∣∣n2,(5)结合(4),(5)，我们显然有
(6)∣∣Xβ0−Xβ^∣∣n2+2λ∣∣β^−β0∣∣1≤16λ2∣S0∣ϕ^2(3,S0).||X\beta^0 -X\hat{\beta}||_n^2+2\lambda||\hat{\beta}-\beta^0||_1\leq\frac{16\lambda^2|S^0|}{\hat{\phi}^2(3,S^0)}.\tag6∣∣Xβ0−Xβ^∣∣n2+2λ∣∣β^−β0∣∣1≤ϕ^2(3,S0)16λ2∣S0∣.(6)
这就是所谓的oracle不等式，注意到如果C-Constant远离0，且λ≍log⁡pn\lambda\asymp\sqrt{\frac{\log p}{n}}λ≍nlogp我们有∣∣Xβ0−Xβ^∣∣n2=Op(∣S0∣log⁡pn).||X\beta^0 -X\hat{\beta}||_n^2=O_p(\frac{|S^0|\log p}{n}).∣∣Xβ0−Xβ^∣∣n2=Op(n∣S0∣logp).

我这里只是给出了对Oracle的一个直观感受，而书中用了一些特别的技巧，也给出了一个更加严格而且一般的定理：

定理2.2(Oracle)

假设∣∣XTϵ/n∣∣∞≤λ0,0≤δ<1||X^T\epsilon/n||_\infty\leq\lambda^0,0\leq\delta<1∣∣XTϵ/n∣∣∞≤λ0,0≤δ<1且λ>λ0\lambda>\lambda_0λ>λ0，令λ−=λ−λ0,λ−=λ+λ0+δλ−,L=λ−(1−δ)λ−,\lambda_{-}=\lambda-\lambda_0,\lambda^-=\lambda+\lambda_0+\delta\lambda_-,L=\frac{\lambda^-}{(1-\delta)\lambda_-},λ−=λ−λ0,λ−=λ+λ0+δλ−,L=(1−δ)λ−λ−,那么我们有
2δλ−∣∣β^−β∣∣1+∣∣Xβ0−Xβ^∣∣n22\delta\lambda_-||\hat{\beta}-\beta||_1+||X\beta^0 -X\hat{\beta}||_n^22δλ−∣∣β^−β∣∣1+∣∣Xβ0−Xβ^∣∣n2≤min⁡β∈Rpmin⁡S⊂{1,...,p}{2δλ−∣∣β−β0∣∣1+∣∣Xβ0−Xβ∣∣n2\leq\min_{\beta\in\mathbb{R}^p}\min_{S\subset\{1,...,p\}}\{2\delta\lambda_-||\beta-\beta^0||_1+||X\beta^0 -X\beta||_n^2≤β∈RpminS⊂{1,...,p}min{2δλ−∣∣β−β0∣∣1+∣∣Xβ0−Xβ∣∣n2+λ−2∣S∣ϕ^2(L,S)+4λ∣∣β−S∣∣1}.+\frac{{\lambda^-}^2|S|}{\hat{\phi}^2(L,S)}+4\lambda||\beta_{-S}||_1\}.+ϕ^2(L,S)λ−2∣S∣+4λ∣∣β−S∣∣1}.如果(β∗,S∗)(\beta^*,S^*)(β∗,S∗)是不等式右边的一个最小元，那么我们就称(β∗,S∗)(\beta^*,S^*)(β∗,S∗)是一个Oracle。显然，取β=β0\beta=\beta^0β=β0，可以得到我们之前得出的上面那个结果。

参考资料

[1] Sara van de geer, Estimation and Testing Under Sparsity, 2016