M.Su CUFE

LARS算法及代码实现

2018-06-29
Algorithm  R  python  C++
     

LARS算法及代码实现

  • 最小角度回归(Least Angle Regression)是lasso的一种高效解法,它不会像前向梯度（Forward Stagewise）算法一样每次走的太小,迭代次数太多.每向前走一步都确保有变量的进入或退出.
  • 本文简单介绍了lars算法的原理,并且给出了python,R语言以及使用RcppArmadillo三个版本的代码

LARS算法原理

lasso的两个等价形式:

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

记
则有

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

  • 求解要求:

  • 求导有：

  • 由互补松弛定理有:
  • 同理

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

注：

  • $\lambda$和t的变化趋势相反，当t逐渐变大的时候，$\lambda$逐渐变小。

  • 记 当$\lambda\geq\lambda_0$时,$\hat{\beta}(\lambda)=0.$

  • 记活跃集为$A(\lambda)={j|\hat{\beta}_j(\lambda)≠0}$

  • 存在$\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,$ 解得

  • 记活跃集为$A(\lambda)={j|\hat{\beta}_j(\lambda)≠0}$

  • 存在$\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)

  • 可以验证，

  • 式(7)减式(6)，
  • 记
    则

活跃集，考虑变量进入

首先定义以下变量

  • 记$A^C_k$为$A_k$的补集
  • 记 $c_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$保持不变，那么，

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

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

画图找交点

  • 若$c_{kj}>0$
  • 若$c_{kj}<0$
  • 综上，

考虑变量退出

  • $\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$

    则
  • $\forall \gamma_j\leq0$
    则
  • 更新公式

LARS的python实现

定义类LARSF,给定.csv文件路径直接读入数据.

属性 描述
x 自变量矩阵
y 因变量
xx $X^TX$
xy $X^TY$
n 样本个数
p 变量个数
b 系数估计矩阵
lam 正则项lambd
函数 作用 输入 返回
forupdate 当有变量进入时更新$X^TX$的cholesky分解 L:下三角矩阵 更新后的下三角矩阵
    xxk:待更新的$X^TX$的行  
    xkxk:待更新的值  
givens 当有变量退出时更新$X^TX$的cholesky分解 L:下三角矩阵 更新后的下三角矩阵
    k:退出变量的位置  
lar 最小角度回归 lam缺省值为-1 lam不指定或小于0时返回整个估计路径
      输入lam时返回对应lambda的估计 
import numpy as np
from scipy import linalg

class LARSF(object):
    def __init__(self,xurl,yurl):
        self.x=np.loadtxt(xurl,delimiter=',',skiprows=0)
        self.y=np.loadtxt(yurl,delimiter=',',skiprows=0)
        self.n=self.x.shape[0]
        self.p=self.x.shape[1]
        self.xx=np.dot(self.x.T,self.x)
        self.xy=np.dot(self.x.T,self.y)
        self.b = []
        self.lam = []
        
        
        
    def forupdate(self,L,xxk,xkxk):
        lk = linalg.solve(L,xxk)
        lkk = np.sqrt(xkxk - sum(lk * lk))
        lk = np.append(lk,lkk)
        zero = np.zeros(L.shape[0])
        LL = np.column_stack((L,zero))
        L = np.row_stack((LL,lk))
        return L
    
    def givens(self,L,k):
        p = L.shape[0]
        L = np.delete(L,k,axis=0)
        mk = k
        while(mk<p-1):
            mx = L[mk,mk:mk+2].copy()
            lmx = np.sqrt(sum(mx*mx))
            L[mk,mk] = lmx
            L[mk,mk+1] = 0
            gives = np.array([[mx[0]/lmx,-mx[1]/lmx],[mx[1]/lmx,mx[0]/lmx]])
            if(mk<p-2):
                L[(mk+1):(p-1), mk:(mk+2)] = np.dot(L[(mk+1):(p-1), mk:(mk+2)],gives)
            mk = mk+1
        L = np.delete(L,p-1,axis=1)
        return L
    
    def lar(self,lam=-1):
        w = np.sqrt(np.diag(self.xx))
    
        w1 = 1/w
        A = np.array([False]*self.p)
        j = np.argmax(w1*np.abs(self.xy))
        A[j] = ~A[j]
        VA = np.array([j])
        nVA = np.arange(self.p)
        nVA = np.delete(nVA,j)
        L = np.array([w[j]])
        lamb = w1[j]*abs(self.xy[j])
        b = np.zeros(self.p,dtype=float)
        reb = b.copy()
        relamb = np.array([lamb])
        while(True):
            
            CC = w1*(self.xy - np.dot(self.xx ,b))
            
            SCC = np.sign(CC)
            SCCA = SCC[VA]
            
            td = linalg.solve(L,w[VA]*SCCA)
            d = linalg.solve(L.T, td)
           
            a = w1[nVA].reshape((-1,1))*np.dot(self.xx[nVA,:][:,VA].reshape((nVA.size,VA.size)),d.reshape((-1,1)))
            
            gam = np.zeros(self.p,dtype=float)
            ww = -b[VA]/d
           
            for i in range(VA.size):
                if(ww[i] > 0 and ww[i] < lamb):
                    gam[VA[i]] = ww[i]
                else:
                    gam[VA[i]] = lamb
            mm=max(gam[VA])+1
           
            if(sum(A)<self.p):
                for i in range(nVA.size):
                    if(a[i]*lamb <= CC[nVA[i]]):
                        gam[nVA[i]] = (lamb - CC[nVA[i]])/(1-a[i])
                    else:
                        gam[nVA[i]] = (lamb + CC[nVA[i]])/(1+a[i])
                    if(gam[nVA[i]] > 0):
                        gam[nVA[i]] = gam[nVA[i]]
                    else: 
                        gam[nVA[i]] = mm         
                                      
            j = np.argmin(gam)
            
            gammin = gam[j]
            
            if(lam >= 0 and (lamb - gammin) <= lam):
                b[VA] = b[VA] + (lamb - lam) * d
                self.b = b.copy
                return
            b[VA] = b[VA] + gammin*d
            lamb = lamb - gammin
            relamb = np.append(relamb,lamb)
            reb = np.row_stack((reb,b))
            
            if(lamb == 0):
                break
            jj = np.array(np.where(VA==j))
            
            if(jj.size==0):
                XTXAJ = self.xx[VA,j] 
                XTXJJ = self.xx[j,j]
                L = self.forupdate(L,XTXAJ,XTXJJ) 
                VA = np.append(VA,j)
                de = np.where(nVA==j)
                nVA = np.delete(nVA,de[0])
            else:
                
                L = self.givens(L,jj[0][0])
                nVA = np.append(nVA,VA[jj[0][0]])
                VA = np.delete(VA,jj[0][0])
                
            A[j] = ~A[j]
            
        self.b = reb.copy   
        self.lam = relamb.copy
        return

LARS的R语言实现

gives <- function(mx, lmx){  
    mc <- mx[1]/lmx                               
    ms <- mx[2]/lmx                                
    matrix(c(mc,ms,-ms,mc),ncol=2)
    }   
  mgives <- function(L,k){     #mgives用于lars中变量退出时xtx的更新
    p <- dim(L)[1]
    if( k>p ) return ("Wrong input of k!")
    Lk <- L[-k,]
    mk <- k
    while( mk < p ){
      mx <- Lk[mk,mk:(mk+1)]
      lmx <- sqrt(sum(mx*mx))
      Lk[mk,mk:(mk+1)] <- c(lmx,0)
      if( mk < p-1 ){
        Lk[(mk+1):(p-1), mk:(mk+1)] <- Lk[(mk+1):(p-1), mk:(mk+1)] %*% gives(mx, lmx)
      }
      mk <- mk + 1
    }
    return(Lk[,-p])
  }
  forupdate <- function(L, xxk, xkxk){   ####forupdate用于lars中变量进入时xtx的更新
    lk <- forwardsolve(L, xxk)
    lkk <- sqrt(xkxk - sum(lk*lk))
    return( as.matrix( rbind( cbind(L,0),c(lk,lkk) ) ) )
  }
lars <- function(XTX, XTY, lam=NULL){           ####lars输入为中心化以后的xtx,xty,lam为选定的lambda
    p <- nrow(XTX)                                    #当lam不为空时，迭代到lambda=lam时停止迭代，并输出估计系数
    w <- sqrt(diag(XTX))                              #当lam为空时，迭代至0，输出list(relamb,reb)
    w1 <- 1/w
    relamb <- NULL
    reb <- NULL
    A <- rep(F,p)
    VA <- NULL
    j <- which.max((1/w)*abs(XTY))
    A[j] <- !A[j]
    VA <- c(VA,j)
    L <- matrix(w[j],1)
    lamb <- w1[j]*abs(XTY[j])
    b <- rep(0,p)
    relamb <- c(relamb, lamb)
    reb <- rbind(reb ,b)
    while(TRUE){
      CC <- w1*(XTY - XTX %*% b)
      SCC <- sign(CC)
      SCCA <- SCC[VA]
      td <- forwardsolve(L,w[VA]*SCCA)
      d <- backsolve(t(L),td)
      a <- w1[-VA]*drop(matrix(XTX[-VA,VA],ncol = sum(A))%*%matrix(d,ncol=1))
      gam <- rep(0,p)
      ww <- -b[VA]/d
      gam[VA] <- ifelse( ww>0 & ww<lamb, ww, lamb)
      mm <- max(gam[VA])+1
      if(sum(A)<p) {
        gam[-VA] <-ifelse(a*lamb<=CC[-VA], (lamb-CC[-VA])/(1-a),(lamb+CC[-VA])/(1+a))
        gam[-VA] <- ifelse(gam[-VA]>0,gam[-VA],mm)                
        
      }
      j <- which.min(gam)
      gammin <- gam[j]
      b[VA] <-b[VA] + gammin * d
      lamb <- lamb - gammin
      relamb <- c(relamb, lamb)
      reb <- rbind(reb, b)
      if((!is.null(lam))&&lamb<=lam){
        step <- length(relamb)
        b <- apply(reb[(step-1):step,],2,function(yy,xx,xout)approx(xx,yy,xout)$y,relamb[(step-1):step],lam)
        return(b)
      }
      if(lamb==0) break
      jj <- which(as.vector(VA)==j)
      if(length(jj)==0){
        XTXAJ <- XTX[VA,j,drop=T] 
        XTXJJ <- XTX[j,j]
        L <- forupdate(L,XTXAJ,XTXJJ) 
        VA <- c(VA,j)
      }else{
        L <- mgives(L,jj) 
    
        VA <- VA[-jj]
      }
      A[j] <- !A[j]
    }
    list(relamb=relamb, reb=reb)                          
  }

LARS的RcppAramdillo实现

#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
using namespace std;

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins("cpp11")]]

void forupdate(mat &L, vec &xxk,double &xkxk)
{
  vec lk = solve(trimatl(L), xxk);
  double lkk = sqrt(xkxk - sum(lk % lk));
  lk.resize(lk.size()+1);
  lk[lk.size()-1] = lkk;
  vec zero(L.n_rows, fill::zeros);
  L = join_cols(join_rows(L,zero), lk.t());
}
mat givens(mat L, uword& k)
{
  int p = L.n_cols;
  int n = L.n_rows;
  
  
  mat Lk = L;
  Lk.shed_row(k);
  
  int mk = k;
  while(mk < p-1){
    vec mx = {Lk(mk,mk),Lk(mk,mk+1)};
    double lmx = sqrt(sum(mx % mx));
    Lk(mk,mk) = lmx;
    Lk(mk,mk+1) = 0;
    mat gives = {mx(0)/lmx, -mx(1)/lmx},{mx(1)/lmx,mx(0)/lmx};
    
    if(mk < p-2)
      Lk.submat(mk+1,mk,p-2,mk+1) = Lk.submat(mk+1,mk,p-2,mk+1) * gives;
    mk += 1;
  }
  Lk.shed_col(p-1);
  return(Lk);
}
// [[Rcpp::export]]
List lar_cpp(const mat& XTX, const vec& XTY, double lam1=-1)
{ 
  int p = XTX.n_rows;
  vec ww=XTX.diag();
  mat w2 = diagmat(ww);
  vec w = sqrt(ww);
  vec w1 = 1/w;
  LogicalVector A(p);
  A.fill(false);
  
  
  uword j = (abs(XTY)%w1).index_max();
  A[j] = !A[j];
  uvec VA ={j};
  uvec nonVA=linspace<uvec>( 0, p-1 ,p ) ;
  nonVA.shed_row(j);
  mat L = {w[j]};
  double lamb = w1[j]*abs(XTY[j]);
  vec b = zeros(p);
  NumericVector relamb ={lamb};
  mat reb = b.t();
  
  while(TRUE){
    vec CC = w1%(XTY - XTX * b);
    vec SCC = sign(CC);
    vec SCCA = SCC(VA);
    vec td = solve(trimatl(L),w(VA)%SCCA);
    vec d = solve(trimatu(L.t()),td);
    vec a = w1.elem(nonVA)%(XTX.submat(nonVA,VA)*d);
    vec gam = zeros(p);
    vec ww = -b(VA)/d;
    for(int i = 0; i <VA.size();i++){
      if(ww[i] > 0 & ww[i] < lamb)
        gam(VA[i]) = ww[i];
      else
        gam(VA[i]) = lamb;
    }
    if(sum(A)<p) {
      for(int i = 0; i < nonVA.size(); i++){
        if(a[i]*lamb <= CC(nonVA[i]))
          gam(nonVA[i]) = (lamb - CC(nonVA[i]))/(1-a[i]);
        else
          gam(nonVA[i]) = (lamb + CC(nonVA[i]))/(1+a[i]);
      }
    }
    
    j = gam.index_min();
    double gammin = gam[j];
    if(lam1>=0 & (lamb-gammin)<=lam1){
      b(VA) = b(VA) + (lamb-lam1)*d;
      NumericVector NVb(b.n_elem);
      for(uword i=0;i<b.n_elem;i++){
        
        NVb(i)=b(i);
        
      }
      return(List::create(Named("b") = NVb));
    }      
    b(VA) = b(VA) + gammin*d;
    lamb = lamb - gammin;
    relamb.push_back(lamb);
    reb = join_cols(reb,b.t());
    
    if(lamb==0) break;
    uvec jj = find(VA==j);
    
    if(jj.n_elem==0){
      uvec uvecj={j};
      vec XTXAJ = XTX.submat(VA,uvecj);
      double XTXJJ = w2(j,j);
      forupdate(L,XTXAJ,XTXJJ);
      VA.resize(VA.n_elem+1);
      VA(VA.n_elem-1)=j;
      uvec elim=find(nonVA==j);
      nonVA.shed_row(elim[0]);
    }else{
      L = givens(L,jj[0]);
      nonVA.resize(nonVA.n_elem+1);
      nonVA(nonVA.n_elem-1)=VA(jj[0]);
      VA.shed_row(jj[0]);
    }
    A[j] = !A[j];
  }

  return( List::create(_["reb"] = reb , _["relamb"] = relamb));
  
}

Similar Posts

上一篇 最优子集回归

下一篇 python学习笔记--生成器

Comments

Content