最小角回归详解

2021-06-29 本文已影响0人 Boye0212

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

1 算法介绍

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

先对 $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$ ；
找出与 $r$ 相关性最大的 $x_{\cdot j}$ ，加入active set；
将 $\beta_j$ 从 $0$ 逐步向LS系数 $x_{\cdot j}'r$ 变动，直到有另一个 $x_{\cdot k}$ ，它与 $r$ 的相关系数绝对值，和 $x_{\cdot j}$ 与 $r$ 的相关系数绝对值一样大；
将 $\beta_j$ 和 $\beta_k$ 同时向二者的联合LS系数变动，直到再出现下一个 $x_{\cdot l}$ ，它与 $r$ 的相关系数满足上一步的条件；
重复上述过程， $\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.