最小角回归详解

本文介绍LAR（Least angle regression，最小角回归），由Efron等（2004）提出。这是一种非常有效的求解LASSO的算法，可以得到LASSO的解的路径。

1 算法介绍

我们直接看最基本的LAR算法，假设有(N)个样本，自变量是(p)维的：

1. 先对(X)（(N imes p)）做标准化处理，使得每个predictor（(X)的每列）满足(x_{cdot j}' 1_N=0)，(Vert x_{cdot j}Vert=1)。我们先假设回归模型中只有截距项，则(eta_0=dfrac{1}{N} y' 1_N)，记残差(r=y-1_N eta_0)，而其他的系数(eta_1=cdots=eta_p=0)；
2. 找出与(r)相关性最大的(x_{cdot j})，加入active set；
3. 将(eta_j)从(0)逐步向LS系数(x_{cdot j}'r)变动，直到有另一个(x_{cdot k})，它与(r)的相关系数绝对值，和(x_{cdot j})与(r)的相关系数绝对值一样大；
4. 将(eta_j)和(eta_k)同时向二者的联合LS系数变动，直到再出现下一个(x_{cdot l})，它与(r)的相关系数满足上一步的条件；
5. 重复上述过程，(min(N-1,p))步后，就得到完整的LS解。

2 算法性质

2.1 保持最小角

我们先来看LS估计量的一个性质：若每个predictor与(y)的相关系的数绝对值相等，从此时开始，将所有系数的估计值同步地从(0)移向LS估计量，在这个过程中，每个predictor与残差向量的相关系数会同比例地减少。

假设我们标准化了每个predictor和(y)，使他们均值为(0)，标准差为(1)。在这里的设定中，对于任意(j=1,ldots,p)，都有(left|x_{cdot j}'y ight|/N=lambda)，其中(lambda)为常数。LS估计量(hateta=(X'X)^{-1}X'y)，当我们将系数从(0)向(hateta)移动了(alpha)（(alphain[0,1])）比例时，记拟合值为(u(alpha)=alpha Xhateta)。

另外，记(ell_p^{(j)})为只有第(j)个元素为(1)、其他元素均为(0)的(p)维向量，则(x_{cdot j}=Xell_p^{(j)})，再记( ext{RSS}=Vert y-XhatetaVert^2)，记投影矩阵(P=X(X'X)^{-1}X')。

这里的问题是，在(alpha)变大过程中，每一个(x_{cdot j})与新的残差的相关系数，是否始终保持相等？且是否会减小？

由于(left| x_{cdot j}' [y-u(alpha)] ight|=left|x_{cdot j}'y - ell_p^{(j)prime} X' u(alpha) ight|=(1-alpha)Nlambda)，即内积与(j)无关。再由( ext{RSS}=(y-Py)'(y-Py)=N-y'Py)可知(y'Py=N- ext{RSS})。

相关系数的绝对值

[egin{aligned} lambda(alpha)=& dfrac{left| x_{cdot j}' [y-u(alpha)] ight|}{Vert x_{cdot j}Vert Vert y-u(alpha)Vert}\ =& dfrac{(1-alpha)Nlambda}{sqrt{N} sqrt{[y-u(alpha)]'[y-u(alpha)]}}\ =& dfrac{(1-alpha)Nlambda}{sqrt{N} sqrt{N(1-alpha)^2+(2alpha-alpha^2) ext{RSS}}}\ =& egin{cases} dfrac{lambda}{sqrt{1+left[-1+dfrac{1}{(1-alpha)^2} ight]dfrac{ ext{RSS}}{N}}},&alphain [0,1)\ 0,&alpha=1 end{cases} end{aligned} ]

因此，任意predictor与当前残差的相关系数绝对值，会随着(alpha)的增加，同比例地减小，并且(lambda(0)=lambda)，(lambda(1)=0)。

现在，我们再回顾一下LAR的过程。在第(k)步开始时，将所有active set中的predictor的集合记为(mathcal{A}_k)，此时在上一步估计完成的系数为(hateta_{mathcal{A}_k})，它是(k-1)维且每个维度都非零的向量，记此时残差为(r_k=y-X_{mathcal{A}_k}hateta_{mathcal{A}_k})，用(r_k)对(X_{mathcal{A}_k})做回归后系数为(delta_k=(X_{mathcal{A}_k}'X_{mathcal{A}_k})^{-1}X_{mathcal{A}_k}' r_k)，拟合值(u_k=X_{mathcal{A}_k}delta_k)。另外，我们知道(X_{mathcal{A}_k}'u_k=X_{mathcal{A}_k}'r_k)，而一个predictor加入(mathcal{A}_k)的条件就是它与当前(r_k)的相关系数的绝对值等于(mathcal{A}_k)中的predictor与当前(r_k)的相关系数的绝对值，所以(X_{mathcal{A}_k}' r_k)向量的每个维度的绝对值都相等，也即(X_{mathcal{A}_k}' u_k)的每个维度的绝对值都相等，(u_k)就是与各个(mathcal{A}_k)中的predictor的角度都相等的向量，且与它们的角度是最小的，而(u_k)也是下一步系数要更新的方向，这也是“最小角回归”名称的由来。

2.2 参数更新

那么，在这个过程中，是否需要每次都逐步小幅增加(alpha)，再检查有没有其他predictor与残差的相关系数绝对值？有没有快速的计算(alpha)的方法？答案是有的。

在第(k)步的开始，(mathcal{A}_k)中有(k-1)个元素，我们记(hat c=X'r_k)，其中(r_k=y-hat y_{mathcal{A}_k})，并记(hat C=max_j {left|hat c_j ight|})，此时的active set其实就是(mathcal{A}_k={j:left|hat c_j ight|=hat C})。在这里，我们将(X_{mathcal{A}_k})做个修改，记(s_j= ext{sign}(hat c_j))，再令(X_{mathcal{A}_k}=[cdots s_jx_{cdot j}cdots]_{jinmathcal{A}_k})。

此时更新方向为(u_k)，(X_{mathcal{A}_k}' u_k=1_{k-1}hat C)，并取(aequiv X' u_k)。更新的规则为(hat y_{mathcal{A}_k}(alpha)= hat y_{mathcal{A}_k}+alpha u_k)。因此，任一predictor，与当前残差的内积就为(c_j(alpha)=hat c_j-alpha a_j)，而对于(jin mathcal{A}_k)，有(left| c_j(alpha) ight|=hat C-alpha hat C)。

对于(jin mathcal{A}_k^c)，如果要使(x_{cdot j})与当前残差的相关系数绝对值，与在(mathcal{A}_k)中的predictor与当前残差的相关系数绝对值相等，也即它们的内积的绝对值相等，必须要满足(|hat c_j-alpha a_j|=(1-hatalpha_j)hat C)。问题转化为了求解使它们相等的(hatalpha_j)，并对于所有的(jin mathcal{A}_k^c)，最小的(hatalpha_j)即为最后的更新步长。

由于(|hat c_j|lt hat C)，因此只需考虑(hat c_j)与(a_j)的大小关系即可。最后解为

[hatalpha_j=egin{cases} dfrac{hat C-hat c_j}{hat C-a_j}, & hat c_jgt a_j\ dfrac{hat C+hat c_j}{hat C+a_j}, & hat c_jleq a_j\ end{cases} ]

注意到

[dfrac{hat C-hat c_j}{hat C-a_j}-dfrac{hat C+hat c_j}{hat C+a_j}=dfrac{2hat C(a_j-hat c_j)}{hat C^2-a_j^2} ]

因此，当(hat c_jgt a_j)时，除非(a_jlt -hat C)即(dfrac{hat C+hat c_j}{hat C+a_j}lt 0)，否则必有(dfrac{hat C-hat c_j}{hat C-a_j} lt dfrac{hat C+hat c_j}{hat C+a_j})。反之，当(hat c_jleq a_j)时，除非(a_jgt hat C)即(dfrac{hat C-hat c_j}{hat C-a_j}lt 0)，否则必有(dfrac{hat C-hat c_j}{hat C-a_j} geq dfrac{hat C+hat c_j}{hat C+a_j})。综上所述，上面的解可以写为

[hat alpha=min_{jin mathcal{A}_k^c}left{dfrac{hat C-hat c_j}{hat C-a_j},dfrac{hat C+hat c_j}{hat C+a_j} ight}^+ ]

其中({}^+)表示只对其中正的元素有效，而丢弃负的元素。

3 LAR与LASSO

LAR虽然是求解LASSO的算法，但它得到的解的路径，在出现了某个系数要穿过(0)的情况时，有可能与LASSO不一样。因此，想要完全得到LASSO的解的路径，还需要做修正。

我们在第1节算法的第4步中加入一个规则：

• 若一个非零系数又变为了(0)，将该predictor从active set中剔除，重新计算当前的LS解作为更新方向。

在修正后，LAR就可以解任意LASSO问题，包括(pgg N)的问题。

为什么会出现与LASSO解不同的情况？我们注意到，对于LASSO的active set (mathcal{B})中的predictor，它的系数需要满足

[x_{cdot j}'(y-Xhateta) = lambda ext{sign}(hateta_j) ]

而对于LAR的active set (mathcal{A})中的predictor，它的系数需要满足

[x_{cdot j}'(y-Xhateta) = gamma s_j ]

其中(s_j)为左边内积的符号。

在正常情况下，上面二者的右侧是相等的，也因此LAR的解就是LASSO的解。但是，当一个非零系数要穿过(0)时，它不再满足LASSO的解条件，因此会被踢出(mathcal{B})，而LAR的解条件却可能没有突变（因为(s_j)是由内积的符号而非系数的符号决定的）。在系数到达(0)时，它满足

[x_{cdot j}'(y-Xhateta) leq lambda ]

这恰恰与(mathcal{A}^c)中的predictor的条件一致，因此可以将它也踢出(mathcal{A})，这样就让LAR与LASSO相一致了。

参考文献

• Efron, Bradley, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. "Least angle regression." Annals of statistics 32, no. 2 (2004): 407-499.
• Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction. Springer Science & Business Media, 2009.