lasso算法.doc

LARS算法简介2011/04/23 郝智恒常规引用方式郝智恒. LARS算法简介. 统计之都, 2011.04. URL: /2011/04/an-introduction-to-lars/.BibTeX引用ARTICLE, AUTHOR = 郝智恒, TITLE = LARS算法简介, JOURNAL = 统计之都, YEAR = 2011, month = 04, URL = /2011/04/an-introduction-to-lars/,最近临时抱佛脚，为了讨论班报告Group Regression方面的文章，研究了Efron等人于2004年发表在Annals of Statistics里一篇被讨论的文章LEAST ANGLE REGRESSION。这篇文章很长，有45页。加上后面一些模型方面大牛的讨论的文章，一共有93页。对于这种超长论文，我向来敬畏。后来因为要报告的文章里很多东西都看不懂，才回过头来研读这篇基石性的文章。所谓大牛，就是他能提出一种别人从来没有提出过的想法。大牛们看待问题的角度和常人不同。比如在回归中常用的逐步回归法。我们小辈们只知道向前回归，向后回归还有二者结合的一些最基本的想法。比如向前回归，就是先选择和响应最相关的变量，进行最小二乘回归。然后在这个模型的基础上，再选择和此时残差相关度最高的（也就是相关度次高）的变量，加入模型重新最小二乘回归。之后再如法继续，直到在某些度量模型的最优性准则之下达到最优，从而选取一个最优的变量子集进行回归分析，得到的模型是相比原模型更加简便，更易于解释的。这种方法，牺牲了模型准确性（预测有偏），但是提高了模型的精确度（方差变小）。大多数本科生对逐步回归的理解也就如此了。Efron看待这个问题时，比起常人更高了一个层次。他首先指出，逐步向前回归，有可能在第二步挑选变量的时候去掉和X1相关的，但是也很重要的解释变量。这是因为它每次找到变量，前进的步伐都太大了，侵略性太强。因此在这个基础上，Efron提出了Forward stagewise。也就是先找出和响应最相关的一个变量，找到第一个变量后不急于做最小二乘回归，而是在变量的solution path上一点一点的前进(所谓solution path是指一个方向，逐步回归是在这个方向上进行)，每前进一点，都要计算一下当前的残差和原有的所有变量的相关系数，找出绝对值最大的相关系数对应的变量。我们可以想像，刚开始，前进的步伐很小，相关系数绝对值最大的对应的变量一定还是第一步选入的变量。但是随着前进的进程不断向前，这个相关系数的绝对值是在慢慢减小的，直到找到另外一个变量X2，它和当前前残差的相关系数和第一个入选变量X1的相关系数绝对值相同，并列第一。此时把X2也加入回归模型中，此时回归模型在X1上的系数已经确定了，如果在X1的solution path上继续前进，则得到的与当前残差相关系数最大的变量一定是X2，所以不再前进，而是改为在X2的solution path上前进，直到找到第三个变量X3，使得X3的与当前残差的相关系数绝对值最大。这样一步一步进行下去。每一步都是很多小步组成。直到某个模型判定准则生效，停止这个步骤。在每一个solution path上的计算都是线性的。总体的solution path是分段线性的。这种算法是一种自动进行模型构建的方法。它和传统的Forward selection在本质上是一样的，都是选择一个变量，然后选择一个继续进行的solution path，在该方向上前进。这两种方法的solution path的选择方法是一样的，唯一的区别就是前进的步伐不一样，Forward selection的前进步伐很大，一次到头，而stagewise则是一小步一小步前进。这样比Forward selection要谨慎一些，会免于漏掉一些重要的变量。从这个视角来看，我们可以选择另外一种solution path。Efron等人在这篇文章中，就提出了一种新的solution path。在已经入选的变量中，寻找一个新的路径，使得在这个路径上前进时，当前残差与已入选变量的相关系数都是相同的。直到找出新的与当前残差相关系数最大的变量。从几何上来看，当前残差在那些已选入回归集的变量们所构成的空间中的投影，是这些变量的角平分线。下面我简单的描述一下这个算法： 第一步，我们初始的估计模型为0，那么当前的残差就是Y,我们找出XY中绝对值最大的那个对应的变量，记为X1,把它加入回归模型。这一步中XY是当前残差和所有变量的相关系数向量。（注意这里Y都已经中心化，X中心标准化过了）。 第二步，在已选的变量的solution path上前进，solution path就是s1*X1，s1是X1与当前残差的相关系数的符号。在这个path上前进，直到另外一个变量出现，使得X1与当前残差的相关系数与它和当前残差的相关系数相同。记这个变量为X2，把它加入回归模型中。 第三步，找到新的solution path。Efron在文章中提出了一种找出满足LARS条件的solution path的解法。solution path需要使得已选入模型变量和当前残差的相关系数均相等。因此这样的路径选择它的方向很显然就是的指向（因为的元素都相同，保证了LARS的要求，当然这里或许会有一些其他的解，也能满足LARS的要求，有没有达人能想到或许证明这个解是唯一的）。只要再标准化这个向量，我们便就找到了solution path的方向。在这个方向上前进，直到下一个满足与当前残差相关系数绝对值最大的变量出现。如此继续下去。 LARS算法，保证了所有入选回归模型的变量在solution path上前进的时候，与当前残差的相关系数都是一样的。这一点，比起Forward stagewise要捷径一些，走得更快一些。LARS算法已经在SAS和R中实现了。作为回归模型选择的一种重要的算法，LARS相比起传统的Forward selection和Forward stagewise，既不那么富于侵略性，又比较走捷径。LARS算法在lasso 估计的求解中也有非常好的应用。在Efron等人的同篇论文中有详细的讨论。关于lasso和它的LARS算法，笔者将在今后的文章中介绍。修正的LARS算法和lasso2011/04/25郝智恒在小弟的上一篇文章中，简单的介绍了LARS算法是怎么回事。主要参考的是Efron等人的经典文章least angle regression。在这篇文章中，还提到了一些有趣的看法，比如如何用LARS算法来求解lasso estimate和forward stagewise estimate。这种看法将我对于模型选择的认识提升了一个层次。在这个更高的层次下看回归的变量选择过程，似乎能有一些更加创新的想法。lasso estimate的提出是Tibshirani在1996年RSSB上的一篇文章Regression shrinkage and selection via lasso。所谓lasso，其全称是least absolute shrinkage and selection operator。其想法可以用如下的最优化问题来表述：在限制了的情况下，求使得残差平方和达到最小的回归系数的估值。我们熟悉如何求解限制条件为等号时，回归方程的求解。也就是用lagrange乘子法求解。但是对于这种，限制条件是不等号的情况，该如何求解，则有两种想法。第一种，也是我比较倾向于的方法，是利用计算机程序，对从开始，不断慢慢增加它的值，然后对每个，求限制条件为等号时候的回归系数的估计，从而可以以的值为横轴，作出一系列的回归系数向量的估计值，这一系列的回归系数的估计值就是lasso estimation。另外一种想法，是借助与最优化问题中的KKT条件，用某个黑箱式的算法，求解。（本人对于最优化方面的东西实在是不很熟悉，故不在此弄斧，只求抛砖引玉，能有高手给出这种想法的具体介绍。）lasso estimate具有shrinkage和selection两种功能，shrinkage这个不用多讲，本科期间学过回归分析的同学应该都知道岭估计会有shrinkage的功效，lasso也同样。关于selection功能，Tibshirani提出，当值小到一定程度的时候，lasso estimate会使得某些回归系数的估值是，这确实是起到了变量选择的作用。当不断增大时，选入回归模型的变量会逐渐增多，当增大到某个值时，所有变量都入选了回归模型，这个时候得到的回归模型的系数是通常意义下的最小二乘估计。从这个角度上来看，lasso也可以看做是一种逐步回归的过程。在我的上一篇文章中，提到了Efron对于逐步回归的一种看法，就是在某个标准之下（比如LARS的标准就是要保证当前残差和已入选变量之间的相关系数相等，也就是当前残差在已入选变量的构成空间中的投影，是那些变量的角平分线）选择一条solution path，在这个solution path上proceed，不断吸收新的变量进入，然后调整solution path 继续proceed。那么对于求解lasso的算法，也有一个相应的对应。Efron提出了一种修正的LARS算法，可以用修正的LARS算法来求解所有的lasso estimates。下面我介绍一下这种修正的LARS算法。首先假设我们已经完成了几步LARS steps。这时候，我们已经有了一个回归变量集，我们记这个回归变量集为。这个集合就对应着一个对于的估计，我们记为。这个估值对应着一个lasso方法对于响应的估值（这里我认为LARS估值和lasso估值应该是一样的），lasso的估值，对应着回归系数的lasso估值，回归系数向量的lasso估值我们记为。为了继续进行下一步，我们先给出一个向量的表达式，然后再解释一下它.就是LARS算法的在当前回归变量集下的solution path。那么我们可以把作为的proceed的path。Efron定义了一个向量,这个向量的元素是,其中是入选变量与当前残差的相关系数的符号，也是的符号。对于没有入选的变量，他们对应在中的元素为0.也就是对应着,我们有。将LARS的solution path对应到lasso estimate的path上，这种对应的想法非常值得借鉴。很显然，会在处变号。那么对于我们已经有的lasso estimate,它中的元素会在最小的的那个大于的处变号。我们记之为。如果没有大于,那么就记为无穷大。对于LARS本身而言，在已经有了如今的回归变量集和当前残差的基础上，我们就会有条solution path，在这个solution path上proceed的最大步记为.通过比较和就会有进一步的想法。Efron的文章证明了如果小于，则对应于LARS估计的那个不会成为一个lasso estimation。（这个是因为当前残差和对应变量的相关系数的符号一定是和该变量的系数符号一致才行）。在这种情况下，我们就不能继续在LARS的solution path上继续前进了，为了利用LARS算法求得lasso estimate，Efron提出把所对应的那个所对应的从回归变量中去掉。去掉之后再计算当前残差和当前这些变量集之间的相关系数，从而确定一条新的solution path，继续进行LARS step。这样进行下去，可以通过LARS算法得到所有的lasso estimate。这个对于LARS的lasso修正算法，被Efron称作“one at a time”条件，也就是每一步都要增加或删掉一个变量。下图显示了用修正了的LARS算法求lasso estimate的过程。这个图是Efron等人的文章中，对于一个实际数据进行回归得到的。该数据一共有10个变量。图的横轴，是所有回归系数估值的绝对值之和，这个值从增加。左侧的纵轴，是回归系数的估值，右侧纵轴是这些回归系数对应的变量的下标。这个图中，我们可以看到每一个回归系数的path。可以看到第七个变量对应的回归系数在横轴快到3000的时候变为了0，说明到这一步时，该变量被删除掉，之后又被重新添加到了回归变量集中。下面通过一个简单的模拟，对lars和lasso以及forward stagewise做一个简单的实现。其实在R中已经有了一个名为lars的包，可以实现上述三种回归。首先，我要模拟的方程为其中和是服从二维联合正态分布，均值为零向量，,服从。我取了50次观测，然后分别通过lasso，lars，以及forward stagewise三种算法进行了回归，其变量的回归路径如下图。简单的代码我直接贴在本文的最后。从这三个算法的图中，我们并