2019-04-16

数据挖掘|LARSN

最后编辑于190514

本文需要一定的Lasso基础，可以推荐的文章学习。

维数灾难

高维数据

何谓高维数据？高维数据指数据的维度很高，甚至远大于样本量的个数。高维数据的明显的表现是：在空间中数据是非常稀疏的，与空间的维数相比样本量总是显得非常少。
在分析高维数据过程中碰到最大的问题就是维数的膨胀，也就是通常所说的“维数灾难”问题。研究表明，随着维数的增长，分析所需的空间样本数会呈指数增长。
如下所示，当数据空间维度由1增加为3，最明显的变化是其所需样本增加；换言之，当样本量确定时，样本密度将会降低，从而样本呈稀疏状态。假设样本量n=12,单个维度宽度为3，那在一维空间下，样本密度为12/3=4，在二维空间下，样本分布空间大小为3×3，则样本密度为12/9=1.33，在三维空间下样本密度为12/27=0.44。

设想一下,当数据空间为更高维时，$X=[x_1x_1,x_2x_2,….,x_nx_n]$会怎么样？

需要更多的样本，样本随着数据维度的增加呈指数型增长；
数据变得更稀疏，导致数据灾难；
在高维数据空间，预测将变得不再容易；
导致模型过拟合。

具体例子可以参考机器学习:分类问题中的“维数灾难”

高维数据对非参数方法带来了致命打击。

理论基础：数据降维

对于高维数据，维数灾难所带来的过拟合问题，其解决思路是：
1）增加样本量；
2）减少样本特征。

而对于现实情况，会存在所能获取到的样本数据量有限的情况，甚至远小于数据维度，即：d>>n。如证券市场交易数据、多媒体图形图像视频数据、航天航空采集数据、生物特征数据等。

主成分分析作为一种数据降维方法，其出发点是通过整合原本的单一变量来得到一组新的综合变量，综合变量所代表的意义丰富且变量间互不相关，综合变量包含了原变量大部分的信息，这些综合变量称为主成分。主成分分析是在保留所有原变量的基础上，通过原变量的线性组合得到主成分，选取少数主成分就可保留原变量的绝大部分信息，这样就可用这几个主成分来代替原变量，从而达到降维的目的。

但是，主成分分析法只适用于数据空间维度小于样本量的情况，当数据空间维度很高时，将不再适用。

Lasso是另一种数据降维方法，该方法不仅适用于线性情况，也适用于非线性情况。Lasso是基于惩罚方法对样本数据进行变量选择，通过对原本的系数进行压缩，将原本很小的系数直接压缩至0，从而将这部分系数所对应的变量视为非显著性变量，将不显著的变量直接舍弃。

lasso回归

普通线性模型

$Y = X\beta + \epsilon$

响应变量$Y = (y_1 + y_2 + … + y_n)^T$,协变量$X = (x^{(1)} + x^{(2)} + … + x^{(n)})$，对于每一个$X^{(j)}$有$X^{(j)} = (x^{(j)}_1 + x^{j}_2 + … + x^{(j)}_n)^T$.

假设每个$X^{(j)} (i = 1,2,…,n, j = 1,2,…,d)$均中心化和标准化（实际操作中这一步很关键），随机误差项$\epsilon_i … N(0,\sigma^2),(i = 1,2,…,n),\epsilon = (\epsilon_1,\epsilon_2,…,\epsilon_n)^T$ 回归系数$\beta = (\beta_1,\beta_2,…,\beta_d)^T$

当X为列满秩设计矩阵时，回归系数$\beta$可由普通最小二乘估计方法求得，

$\hat{\beta} = (X^TX)^{-1}X^TY$

当X不为列满秩设计矩阵时，普通最小二乘法将不再适用于求解回归系数$\beta$，惩罚方法应用较广。

惩罚方法

当X不为列满秩设计矩阵时，普通最小二乘法将不再适用于求解回归系数，此时引入惩罚方法。该方法可以同时实现变量选择和参数估计，在参数估计时，通过将部分参数压缩为0来达到变量选择的效果。惩罚方法时取惩罚似然函数最小时的值作为回归系数的估计值，即 $\hat{\beta} = arg min_{\beta∈R^d}(||Y - X\beta||^2 + P_{\lambda}(|\beta|))$

其中惩罚项$P_{\lambda}(|\beta|) = \lambda \sum_{j = 1}^{d}|\beta_j|^m,m>=0$,$\lambda$为调节参数（也可以为向量）。m = 1时，得到$L_1$惩罚项（Lasso惩罚）；当m = 2时，得到$L_2$惩罚项（Ridge惩罚）。

Lasso方法

Lasso方法是在普通线性模型中增加$L_1$惩罚项，对于普通线性模型的Lasso估计为 $\hat{\beta_{Lasso}} = arg min_{\beta∈R^d}||Y - X\beta||^2 s.t. \sum_{j = 1}^{d}|\beta_j|<=t,t>=0$

等价于 $\hat{\beta_{Lasso}} = arg min_{\beta∈R^d}||Y - X\beta||^2 +\lambda \sum_{j = 1}^{d}|\beta_j|$

其中，t与$\lambda$一一对应，为调节系数。

令 $t_0 = \sum_{j = 1}^{d}|\hat{\beta_j}(OLS)|$ ，当t<t0时，一部分系数就会被压缩至0，从而降低X的维度，达到减小模型复杂度的目的。

参考文章Lasso回归

LARS

一般性质

lasso的两个等价形式:

$\hat{\beta}(t) = arg\ min \frac 12\|Y-X\beta||_2^2\ \ \ \ s.t \ \ ||\beta||_1\leq t \tag{1}$ $\hat{\beta}(\lambda) = arg\ min \frac 12 \ ||Y-X\beta||_2^2+\lambda\ ||\beta||_1 \tag{2}$

$\forall t \geq 0,\exists \lambda \geq 0,$ 使得$ \hat{\beta}(\lambda)=\hat{\beta}(t).$

记

$\beta_j^+=\left\{\begin{array}{ll}\beta_j,& \beta_j\geq0\\ 0,& \beta_j< 0\end{array}\right., \quad\quad \beta_j^-=\left\{\begin{array}{ll}-\beta_j,& \beta_j\leq0\\ 0,& \beta_j>0\end{array}\right.$

则有

$\beta_j=\beta_j^+-\beta_j^- ,\quad |\beta_j|=\beta_j^++\beta_j^-$

在考虑引入正部和负部后，增加了两个约束条件，故拉格朗日函数为：

$\Gamma(\beta_j^+,\beta_j^-;\lambda)=\frac 12||Y-\sum(X_j(\beta_j^+-\beta_j^-))||_2^2+\lambda1_p^T(\beta_j^++\beta_j^-)-\sum{\lambda_j^+\beta_j^+}-\sum{\lambda_j^-\beta_j^-}\tag{3}$

求解要求: $\begin{cases}\lambda_j^+\beta_j^+=0\\\lambda_j^-\beta_j^-=0\\ \beta_j^+\geq0 \\ \beta_j^-\geq0 \end{cases}$
求导有： $\begin{cases}\frac{\partial \Gamma(\beta_j^+,\beta_j^-;\lambda)}{\partial \beta_j^+}=-X_j^T(Y-X\beta)+\lambda-\lambda_j^+=0 \\ \frac{\partial \Gamma(\beta_j^+,\beta_j^-;\lambda)}{\partial \beta_j^-}=X_j^T(Y-X\beta)+\lambda-\lambda_j^-=0\end{cases} \tag{4}$
由互补松弛定理有:
$\begin{cases}\lambda_j^+= 0,\quad \beta_j^+>0\\ \lambda_j^+\geq0,\quad \beta_j^+=0 \end{cases}$
同理 $\begin{cases}\lambda_j^-= 0,\quad \beta_j^->0\\ \lambda_j^-\geq0,\quad \beta_j^-=0 \end{cases}$

令导数为0，于是式(4)退化为：

$\begin{cases}若\beta_j^+(\lambda)\neq 0 ,\ \ \ X_j^T(Y-X\beta(\lambda))=\lambda \\ 若\beta_j^-(\lambda)\neq 0 ,\ \ \ X_j^T(Y-X\beta(\lambda))=-\lambda \\ 若\beta_j(\lambda)=0,\ \ \ |X_j^T(Y-X\beta(\lambda))|\leq\lambda\end{cases}\tag{5}$

注：

$\lambda$和t的变化趋势相反，当t逐渐变大的时候，$\lambda$逐渐变小。
记
$\lambda_0= \max\{|X_1^TY|,\cdots,|X_p^TY|\}, \quad j_0 = \arg\max\{|X_1^TY|,\cdots,|X_p^TY|\},$
当$\lambda\geq\lambda_0$时,$\hat{\beta}(\lambda)=0.$
记活跃集为$A(\lambda)=\{j\ |\hat{\beta}_j(\lambda)\neq0\}$
存在$\lambda_k$, $k=0,\dots,$ 当$\lambda\in(\lambda_{k+1},\lambda_{k}),$
$A(\lambda)$保持不变
当$\lambda =0$，该模型的解为最小二乘解

初始化

$\exists\lambda_0$，当$\lambda\geq\lambda_0$时,$\hat{\beta}(\lambda)=0,$
解得 $\lambda_0= \max\{|X_1^TY|,\cdots,|X_p^TY|\}, \quad j_0 = \arg\max\{|X_1^TY|,\cdots,|X_p^TY|\}.$

记活跃集为$A(\lambda)=\{j\ |\hat{\beta}_j(\lambda)\neq0\}$
存在$\lambda_k$, $k=0,\dots,$ 当$\lambda\in(\lambda_{k+1},\lambda_{k}),$
$A(\lambda)$保持不变
在$\lambda=\lambda_k$时，活跃集发生变化，或者有变量进入，或者有变量退出，但可以证明以下事实
- $\hat{\beta}(\lambda)$是连续的
- 当$\lambda\in(\lambda_{k+1},\lambda_{k}),$ $\hat{\beta}(\lambda)$是分段线性的
所以，最小角度回归算法有两个关键点
- 当$\lambda\in(\lambda_{k+1},\lambda_{k})$时，其线性形式是什么
- 确定$\lambda_{k+1},$ 也就是分段线性形式改变的位置

一般化、第k步

当$\lambda\in(\lambda_{k+1},\lambda_{k}),$ $A(\lambda)$保持不变, 记做$A_k$
记$\beta_{A_k}=\beta_{A_k}(\lambda_k),$ $c_{k}=X_{A_K}^T(Y-X_{A_K}\beta_{A_K}), $
$s_{A_K}=sign(c_{k})$
由式(5) $X_{A_K}^T(Y-X_{A_K}\hat{\beta}_{A_K}(\lambda))=s_{A_K}\lambda, \quad \quad\lambda\in(\lambda_{k+1},\lambda_k) \tag{6}$
可以验证，
$X_{A_K}^T(Y-X_{A_K}\hat\beta_{A_K})=s_{A_K}\lambda_k\tag{7}$
式(7)减式(6)，
$\hat\beta_{A_k}(\lambda)=\hat\beta_{A_k}+(\lambda_k-\lambda)(X_{A_k}^TX_{A_k})^{-1}S_{AK}$
记 $d_{A_k}=(X_{A_k}^TX_{A_k})^{-1}s_{A_k}\quad\quad \gamma=\lambda_k-\lambda$
则 $\hat\beta_{A_k}(\lambda)=\hat\beta_{A_k}+\gamma d_{A_k}$

活跃集，考虑变量进入

首先定义以下变量

记$A^C_k$为$A_k$的补集
记 $a_k=(X_{A^C_k}^TX_{A_k})d_{A_k},$ $c_k$为处于非活跃集自变量与残差的相关系数， $c_k=X_{A^C_k}^T(Y-X_{A_k}\hat{\beta}_{A_k})$
记$a_{kj}$是$a_k$的第$j$项, $c_j$是$c$的第j项

那么，

$a_{kj}\lambda_k\leq c_{kj}$时，$\gamma_j=\frac{\lambda_k-c_{kj}}{1-a_{kj}}$
$a_{kj}\lambda_k> c_{kj}$时，$\gamma_j=\frac{\lambda_k+c_{kj}}{1+a_{kj}}$

具体推导如下：

根据KKT条件，如果$A_k$保持不变，那么，

$\left|X_{A^C_k}^T(Y-X_{A_k}\hat\beta_{A_k}({\gamma}))\right|\leq \lambda_k-\gamma,$

式(3)中，$\hat\beta_{A_k}({\gamma})=\hat\beta_{A_k}+\gamma d_{A_k}，$ 所以

$|c_k - \gamma a_k|\leq \lambda_k-\gamma$

如果$A^C_k$中第$j$个元素率先进入活跃集，那么，$c_{kj}-\gamma a_{kj}=\pm(\lambda_k-\gamma)$

画图找交点

若$c_{kj}>0$ $\gamma_j=\begin{cases}\frac{\lambda_k-c_{kj}}{1-a_{kj}},&-a_{kj}\geq -\frac{c_{kj}}{\lambda_k}\\ \frac{\lambda_k+c_{kj}}{1+a_{kj}},&-a_{kj}\leq -\frac{c_{kj}}{\lambda_k}\end{cases}$
若$c_{kj}<0$ $\gamma_j=\begin{cases}\frac{\lambda_k-c_{kj}}{1-a_{kj}},&-a_{kj}\geq -\frac{c_{kj}}{\lambda_k}\\ \frac{\lambda_k+c_{kj}}{1+a_{kj}},&-a_{kj}\leq -\frac{c_{kj}}{\lambda_k}\end{cases}$
综上， $\gamma_j=\begin{cases}\frac{\lambda_k-c_{kj}}{1-a_{kj}},&a_{kj}\lambda_k \leq c_{kj}\\ \frac{\lambda_k+c_{kj}}{1+a_{kj}},&a_{kj}\lambda_k> c_{kj} \end{cases}$

考虑变量退出

$\hat\beta_{A_k}({\gamma})=\hat\beta_{A_k}+\gamma d_{A_k}$，记$w_j$为$\beta_{A_k}$的第$j$个元素与$d_{A_k}$的第$j$个元素之比的相反数
活跃集中第$j$个元素率先退出，那么$\gamma_j= w_j$

活跃集个数小于$p$时，变量的进入和退出

定义$\hat{\gamma}=\min^{+} \{\gamma_j,j=1,2,\cdots,p\}.$

活跃集=p时

$\exists \gamma_j>0$ $记\hat{\gamma}=min{\{\gamma_j,j=1,2,\cdots,p\},其中\gamma_j>0}$ 则 $\hat{\gamma}=\begin{cases}\hat{\gamma},&\lambda_k-\hat{\gamma}>0\\ \lambda_k,&\lambda_k-\hat{\gamma}\leq0 \end{cases}$
$\forall \gamma_j\leq0$
则 $\hat{\gamma}=\lambda_k$
更新公式 $\beta_{new}=\beta_{old}+\hat{\gamma}\cdot d_A$ $\lambda_{new}=\lambda_{old}-\hat{\gamma}$

RSS的更新

$\begin{eqnarray} RSS&=&(Y-X_A(\beta_A+\gamma\cdot d_A))^T(Y-X_A(\beta_A+\gamma\cdot d_A)\\ &=& (Y-X_A\beta_A)^T(Y-X_A\beta_A)-2\gamma(Y-X_A\beta_A)^TX_Ad_A+\gamma^2d_A^TX_A^TX_Ad_A\\ &=& (Y-X_A\beta_A)^T(Y-X_A\beta_A)-2\gamma Y^TX_Ad_A+2\gamma\beta_A^TX_A^TX_Ad_A+\gamma^2d_A^TX_A^TX_Ad_A \end{eqnarray}$

代码实现

l2norm <- function(x){
    return(sqrt(sum(x*x)))
}
l1norm <- function(x){
    return(sum(abs(x)))
}
l0norm <- function(x){
    return( sum(abs(x)>1e-13))
}
    

library(ElemStatLearn)
data(prostate)
data <- prostate[,-10]
head(data)
np <- dim(data)
n <- np[1]
p <- np[2]-1
xn <- names(data)[1:p]
x <- as.matrix(data[,1:p])
xm <- apply(x,2,mean)
X <- sweep(x,2,xm,"-")

Xl2norm <- apply(X,2,l2norm)
X <- sweep(X,2,Xl2norm,"/")
#apply(X,2,sd)
#x <- as.matrix( cbind(1,x) )
y <- data[,p+1]
ym <- mean(y)
Y <- (y-ym)
Yl2norm <- l2norm(Y)
Y <- Y/Yl2norm


##############################
XTX <- t(X)%*%X
XTY <- drop(t(X)%*% Y) 
YY <- sum(Y*Y)
reA <- NULL
relamb <- NULL
reRSS <- NULL
reb <- NULL
########################
A <- rep(F,8)
reA <- rbind(reA, A)
#df <- sum(A)
j <- which.max( abs(XTY))
#print(j)
A[j] <- !A[j]
lamb <- abs(XTY[j])
b <- rep(0,8)
RSS <- YY

relamb <- c(relamb, lamb)
reb <- rbind(reb ,b)
reRSS <- c(reRSS,RSS)
#########################################

############################
while(TRUE){
    CC <- XTY - XTX %*% b
    SCC = sign(CC)
    SCCA <- SCC[A]
    XTXA <- XTX[A,A,drop=F]
    d = as.matrix(solve(XTXA, SCCA),ncol =1)
    a = XTX[!A,A] %*% d
  
    gam=rep(0,8)
    gam[A] = -b[A]/d
    if(sum(A)<=0 | sum(A)>8 ){
        print("Something is wrong!")
        break
    }
    else if(sum(A)<8){
        gam[!A] =ifelse(a*lamb<=CC[!A], (lamb-CC[!A])/(1-a),(lamb+CC[!A])/(1+a))
        gamm= max(gam) +1
        gam[gam<=0]=gamm
  
        j = which.min(gam)
        gammin = gam[j]
        RSS <- RSS - 2*gammin*sum(XTY[A]*d) + 2*gammin*sum(b[A]*SCCA)+
               gammin^2*sum(d*SCCA)
        b[A] =b[A] + gammin * d
        lamb = lamb - gammin
        reA = rbind(reA, A)
        A[j] = !A[j]
        reRSS <- c(reRSS, RSS)
        
        relamb = c(relamb, lamb)
        reb =rbind(reb ,b)
  
    } else if (sum(A)==8) {
        if(sum(gam > 0) > 1) {
            gamm= max(gam) +1
            gam[gam<=0]=gamm
    
            j = which.min(gam)
            gammin = gam[j]
            if(lamb-gammin >0){
                RSS <- RSS - 2*gammin*sum(XTY[A]*d) + 2*gammin*sum(b[A]*SCCA)+
                       gammin^2*sum(d*SCCA)
                b[A] =b[A] + gammin * d
                lamb = lamb - gammin
                reA = rbind(reA, A)
                A[j] = !A[j]
                reRSS <- c(reRSS, RSS)
                reA = rbind(reA, A)
                relamb = c(relamb, lamb)
                reb =rbind(reb ,b)                 
            }else{
                gammin = lamb
                RSS <- RSS - 2*gammin*sum(XTY[A]*d) + 2*gammin*sum(b[A]*SCCA)+
                       gammin^2*sum(d*SCCA)
                
                b[A] =b[A] + gammin * d
                lamb = lamb - gammin
                reRSS <- c(reRSS, RSS)
                reA = rbind(reA, A)
                relamb = c(relamb, lamb)
                reb =rbind(reb ,b)
                break                
            }            
        } else {
             gammin = lamb
             RSS <- RSS - 2*gammin*sum(XTY[A]*d) + 2*gammin*sum(b[A]*SCCA)+
                    gammin^2*sum(d*SCCA)
             
             b[A] =b[A] + gammin * d
             lamb = lamb - gammin
             reRSS <- c(reRSS, RSS)
             reA = rbind(reA, A)
             relamb = c(relamb, lamb)
             reb =rbind(reb ,b)
             break
        }
  }

}
Cp <- reRSS/RSS*(n-p) - n + 2*apply(reA,1,sum)
re <- list(b=reb,lambda=relamb,A=reA,RSS=reRSS,Cp=Cp)
re

结果不贴了，大家跑一下就行啦。

本文链接： https://konelane.github.io/2019/04/16/190416LARSN/

-- EOF --

转载请注明出处署名-非商业性使用-禁止演绎 3.0 国际（CC BY-NC-ND 3.0）

￥^￥请氦核牛饮一盒奶~suki

打赏