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

1 算法介绍

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

  1. 先对\(X\)（\(N\times p\)）做标准化处理，使得每个predictor（\(X\)的每列）满足\(x_{\cdot j}' 1_N=0\)，\(\Vert x_{\cdot j}\Vert=1\)。我们先假设回归模型中只有截距项，则\(\beta_0=\dfrac{1}{N} y' 1_N\)，记残差\(r=y-1_N \beta_0\)，而其他的系数\(\beta_1=\cdots=\beta_p=0\)；
  2. 找出与\(r\)相关性最大的\(x_{\cdot j}\)，加入active set；
  3. 将\(\beta_j\)从\(0\)逐步向LS系数\(x_{\cdot j}'r\)变动，直到有另一个\(x_{\cdot k}\)，它与\(r\)的相关系数绝对值，和\(x_{\cdot j}\)与\(r\)的相关系数绝对值一样大；
  4. 将\(\beta_j\)和\(\beta_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\right|/N=\lambda\)，其中\(\lambda\)为常数。LS估计量\(\hat\beta=(X'X)^{-1}X'y\)，当我们将系数从\(0\)向\(\hat\beta\)移动了\(\alpha\)（\(\alpha\in[0,1]\)）比例时，记拟合值为\(u(\alpha)=\alpha X\hat\beta\)。

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

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

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

相关系数的绝对值

\[\begin{aligned} \lambda(\alpha)=& \dfrac{\left| x_{\cdot j}' [y-u(\alpha)]\right|}{\Vert x_{\cdot j}\Vert \Vert y-u(\alpha)\Vert}\\ =& \dfrac{(1-\alpha)N\lambda}{\sqrt{N} \sqrt{[y-u(\alpha)]'[y-u(\alpha)]}}\\ =& \dfrac{(1-\alpha)N\lambda}{\sqrt{N} \sqrt{N(1-\alpha)^2+(2\alpha-\alpha^2)\text{RSS}}}\\ =& \begin{cases} \dfrac{\lambda}{\sqrt{1+\left[-1+\dfrac{1}{(1-\alpha)^2}\right]\dfrac{\text{RSS}}{N}}},&\alpha\in [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\)，此时在上一步估计完成的系数为\(\hat\beta_{\mathcal{A}_k}\)，它是\(k-1\)维且每个维度都非零的向量，记此时残差为\(r_k=y-X_{\mathcal{A}_k}\hat\beta_{\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\right|\}\)，此时的active set其实就是\(\mathcal{A}_k=\{j:\left|\hat c_j\right|=\hat C\}\)。在这里，我们将\(X_{\mathcal{A}_k}\)做个修改，记\(s_j=\text{sign}(\hat c_j)\)，再令\(X_{\mathcal{A}_k}=[\cdots s_jx_{\cdot j}\cdots]_{j\in\mathcal{A}_k}\)。

此时更新方向为\(u_k\)，\(X_{\mathcal{A}_k}' u_k=1_{k-1}\hat C\)，并取\(a\equiv 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\)，而对于\(j\in \mathcal{A}_k\)，有\(\left| c_j(\alpha)\right|=\hat C-\alpha \hat C\)。

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

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

\[\hat\alpha_j=\begin{cases} \dfrac{\hat C-\hat c_j}{\hat C-a_j}, & \hat c_j\gt a_j\\ \dfrac{\hat C+\hat c_j}{\hat C+a_j}, & \hat c_j\leq 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{2\hat C(a_j-\hat c_j)}{\hat C^2-a_j^2} \]

因此，当\(\hat c_j\gt a_j\)时，除非\(a_j\lt -\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_j\leq a_j\)时，除非\(a_j\gt \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_{j\in \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}\right\}^+ \]

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

3 LAR与LASSO

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

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

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

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

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

\[x_{\cdot j}'(y-X\hat\beta) = \lambda \text{sign}(\hat\beta_j) \]

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

\[x_{\cdot j}'(y-X\hat\beta) = \gamma s_j \]

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

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

\[x_{\cdot j}'(y-X\hat\beta) \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.