本文介绍LAR（Least angle regression，最小角回归），由Efron等（2004）提出。这是一种非常有效的求解LASSO的算法，可以得到LASSO的解的路径。
1 算法介绍
我们直接看最基本的LAR算法，假设有N N N p p p
先对X X X N × p N\times p N × p X X X x ⋅ j ′ 1 N = 0 x_{\cdot j}' 1_N=0 x ⋅ j ′ 1 N = 0 ∥ x ⋅ j ∥ = 1 \Vert x_{\cdot j}\Vert=1 ∥ x ⋅ j ∥ = 1 β 0 = 1 N y ′ 1 N \beta_0=\dfrac{1}{N} y' 1_N β 0 = N 1 y ′ 1 N r = y − 1 N β 0 r=y-1_N \beta_0 r = y − 1 N β 0 β 1 = ⋯ = β p = 0 \beta_1=\cdots=\beta_p=0 β 1 = ⋯ = β p = 0
找出与r r r x ⋅ j x_{\cdot j} x ⋅ j
将β j \beta_j β j 0 0 0 x ⋅ j ′ r x_{\cdot j}'r x ⋅ j ′ r 直到 有另一个x ⋅ k x_{\cdot k} x ⋅ k r r r x ⋅ j x_{\cdot j} x ⋅ j r r r
将β j \beta_j β j β k \beta_k β k 同时 向二者的联合LS系数变动，直到再出现下一个x ⋅ l x_{\cdot l} x ⋅ l r r r
重复上述过程，min ( N − 1 , p ) \min(N-1,p) min ( N − 1 , p )
2 算法性质
2.1 保持最小角
我们先来看LS估计量的一个性质：若每个predictor与y y y 0 0 0
假设我们标准化了每个predictor和y y y 0 0 0 1 1 1 j = 1 , … , p j=1,\ldots,p j = 1 , … , p ∣ x ⋅ j ′ y ∣ / N = λ \left|x_{\cdot j}'y\right|/N=\lambda ∣ ∣ x ⋅ j ′ y ∣ ∣ / N = λ λ \lambda λ β ^ = ( X ′ X ) − 1 X ′ y \hat\beta=(X'X)^{-1}X'y β ^ = ( X ′ X ) − 1 X ′ y 0 0 0 β ^ \hat\beta β ^ α \alpha α α ∈ [ 0 , 1 ] \alpha\in[0,1] α ∈ [ 0 , 1 ] u ( α ) = α X β ^ u(\alpha)=\alpha X\hat\beta u ( α ) = α X β ^
另外，记ℓ p ( j ) \ell_p^{(j)} ℓ p ( j ) j j j 1 1 1 0 0 0 p p p x ⋅ j = X ℓ p ( j ) x_{\cdot j}=X\ell_p^{(j)} x ⋅ j = X ℓ p ( j ) RSS = ∥ y − X β ^ ∥ 2 \text{RSS}=\Vert y-X\hat\beta\Vert^2 RSS = ∥ y − X β ^ ∥ 2 P = X ( X ′ X ) − 1 X ′ P=X(X'X)^{-1}X' P = X ( X ′ X ) − 1 X ′
这里的问题是，在α \alpha α x ⋅ j x_{\cdot j} x ⋅ j
由于∣ x ⋅ j ′ [ y − u ( α ) ] ∣ = ∣ x ⋅ j ′ y − ℓ p ( j ) ′ X ′ u ( α ) ∣ = ( 1 − α ) N λ \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 ∣ ∣ x ⋅ j ′ [ y − u ( α )] ∣ ∣ = ∣ ∣ x ⋅ j ′ y − ℓ p ( j ) ′ X ′ u ( α ) ∣ ∣ = ( 1 − α ) N λ j j j RSS = ( y − P y ) ′ ( y − P y ) = N − y ′ P y \text{RSS}=(y-Py)'(y-Py)=N-y'Py RSS = ( y − P y ) ′ ( y − P y ) = N − y ′ P y y ′ P y = N − RSS y'Py=N-\text{RSS} y ′ P y = N − RSS
相关系数的绝对值
λ ( α ) = ∣ x ⋅ j ′ [ y − u ( α ) ] ∣ ∥ x ⋅ j ∥ ∥ y − u ( α ) ∥ = ( 1 − α ) N λ N [ y − u ( α ) ] ′ [ y − u ( α ) ] = ( 1 − α ) N λ N N ( 1 − α ) 2 + ( 2 α − α 2 ) RSS = { λ 1 + [ − 1 + 1 ( 1 − α ) 2 ] RSS N , α ∈ [ 0 , 1 ) 0 , α = 1 \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} λ ( α ) = = = = ∥ x ⋅ j ∥∥ y − u ( α ) ∥ ∣ ∣ x ⋅ j ′ [ y − u ( α )] ∣ ∣ N [ y − u ( α ) ] ′ [ y − u ( α )] ( 1 − α ) N λ N N ( 1 − α ) 2 + ( 2 α − α 2 ) RSS ( 1 − α ) N λ ⎩ ⎨ ⎧ 1 + [ − 1 + ( 1 − α ) 2 1 ] N RSS λ , 0 , α ∈ [ 0 , 1 ) α = 1 因此，任意predictor与当前残差的相关系数绝对值，会随着α \alpha α λ ( 0 ) = λ \lambda(0)=\lambda λ ( 0 ) = λ λ ( 1 ) = 0 \lambda(1)=0 λ ( 1 ) = 0
现在，我们再回顾一下LAR的过程。在第k k k A k \mathcal{A}_k A k β ^ A k \hat\beta_{\mathcal{A}_k} β ^ A k k − 1 k-1 k − 1 r k = y − X A k β ^ A k r_k=y-X_{\mathcal{A}_k}\hat\beta_{\mathcal{A}_k} r k = y − X A k β ^ A k r k r_k r k X A k X_{\mathcal{A}_k} X A k δ k = ( X A k ′ X A k ) − 1 X A k ′ r k \delta_k=(X_{\mathcal{A}_k}'X_{\mathcal{A}_k})^{-1}X_{\mathcal{A}_k}' r_k δ k = ( X A k ′ X A k ) − 1 X A k ′ r k u k = X A k δ k u_k=X_{\mathcal{A}_k}\delta_k u k = X A k δ k X A k ′ u k = X A k ′ r k X_{\mathcal{A}_k}'u_k=X_{\mathcal{A}_k}'r_k X A k ′ u k = X A k ′ r k A k \mathcal{A}_k A k r k r_k r k A k \mathcal{A}_k A k r k r_k r k X A k ′ r k X_{\mathcal{A}_k}' r_k X A k ′ r k X A k ′ u k X_{\mathcal{A}_k}' u_k X A k ′ u k u k u_k u k A k \mathcal{A}_k A k u k u_k u k
2.2 参数更新
那么，在这个过程中，是否需要每次都逐步小幅增加α \alpha α α \alpha α
在第k k k A k \mathcal{A}_k A k k − 1 k-1 k − 1 c ^ = X ′ r k \hat c=X'r_k c ^ = X ′ r k r k = y − y ^ A k r_k=y-\hat y_{\mathcal{A}_k} r k = y − y ^ A k C ^ = max j { ∣ c ^ j ∣ } \hat C=\max_j \{\left|\hat c_j\right|\} C ^ = max j { ∣ c ^ j ∣ } A k = { j : ∣ c ^ j ∣ = C ^ } \mathcal{A}_k=\{j:\left|\hat c_j\right|=\hat C\} A k = { j : ∣ c ^ j ∣ = C ^ } X A k X_{\mathcal{A}_k} X A k s j = sign ( c ^ j ) s_j=\text{sign}(\hat c_j) s j = sign ( c ^ j ) X A k = [ ⋯ s j x ⋅ j ⋯ ] j ∈ A k X_{\mathcal{A}_k}=[\cdots s_jx_{\cdot j}\cdots]_{j\in\mathcal{A}_k} X A k = [ ⋯ s j x ⋅ j ⋯ ] j ∈ A k
此时更新方向为u k u_k u k X A k ′ u k = 1 k − 1 C ^ X_{\mathcal{A}_k}' u_k=1_{k-1}\hat C X A k ′ u k = 1 k − 1 C ^ a ≡ X ′ u k a\equiv X' u_k a ≡ X ′ u k y ^ A k ( α ) = y ^ A k + α u k \hat y_{\mathcal{A}_k}(\alpha)= \hat y_{\mathcal{A}_k}+\alpha u_k y ^ A k ( α ) = y ^ A k + α u k c j ( α ) = c ^ j − α a j c_j(\alpha)=\hat c_j-\alpha a_j c j ( α ) = c ^ j − α a j j ∈ A k j\in \mathcal{A}_k j ∈ A k ∣ c j ( α ) ∣ = C ^ − α C ^ \left| c_j(\alpha)\right|=\hat C-\alpha \hat C ∣ c j ( α ) ∣ = C ^ − α C ^
对于j ∈ A k c j\in \mathcal{A}_k^c j ∈ A k c x ⋅ j x_{\cdot j} x ⋅ j A k \mathcal{A}_k A k ∣ c ^ j − α a j ∣ = ( 1 − α ^ j ) C ^ |\hat c_j-\alpha a_j|=(1-\hat\alpha_j)\hat C ∣ c ^ j − α a j ∣ = ( 1 − α ^ j ) C ^ α ^ j \hat\alpha_j α ^ j j ∈ A k c j\in \mathcal{A}_k^c j ∈ A k c α ^ j \hat\alpha_j α ^ j
由于∣ c ^ j ∣ < C ^ |\hat c_j|\lt \hat C ∣ c ^ j ∣ < C ^ c ^ j \hat c_j c ^ j a j a_j a j
α ^ j = { C ^ − c ^ j C ^ − a j , c ^ j > a j C ^ + c ^ j C ^ + a j , 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} α ^ j = ⎩ ⎨ ⎧ C ^ − a j C ^ − c ^ j , C ^ + a j C ^ + c ^ j , c ^ j > a j c ^ j ≤ a j 注意到
C ^ − c ^ j C ^ − a j − C ^ + c ^ j C ^ + a j = 2 C ^ ( a j − c ^ j ) C ^ 2 − a j 2 \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} C ^ − a j C ^ − c ^ j − C ^ + a j C ^ + c ^ j = C ^ 2 − a j 2 2 C ^ ( a j − c ^ j ) c ^ j > a j \hat c_j\gt a_j c ^ j > a j a j < − C ^ a_j\lt -\hat C a j < − C ^ C ^ + c ^ j C ^ + a j < 0 \dfrac{\hat C+\hat c_j}{\hat C+a_j}\lt 0 C ^ + a j C ^ + c ^ j < 0 C ^ − c ^ j C ^ − a j < C ^ + c ^ j C ^ + a j \dfrac{\hat C-\hat c_j}{\hat C-a_j} \lt \dfrac{\hat C+\hat c_j}{\hat C+a_j} C ^ − a j C ^ − c ^ j < C ^ + a j C ^ + c ^ j c ^ j ≤ a j \hat c_j\leq a_j c ^ j ≤ a j a j > C ^ a_j\gt \hat C a j > C ^ C ^ − c ^ j C ^ − a j < 0 \dfrac{\hat C-\hat c_j}{\hat C-a_j}\lt 0 C ^ − a j C ^ − c ^ j < 0 C ^ − c ^ j C ^ − a j ≥ C ^ + c ^ j C ^ + a j \dfrac{\hat C-\hat c_j}{\hat C-a_j} \geq \dfrac{\hat C+\hat c_j}{\hat C+a_j} C ^ − a j C ^ − c ^ j ≥ C ^ + a j C ^ + c ^ j
α ^ = min j ∈ A k c { C ^ − c ^ j C ^ − a j , C ^ + c ^ j 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\}^+ α ^ = j ∈ A k c min { C ^ − a j C ^ − c ^ j , C ^ + a j C ^ + c ^ j } + 其中{ } + \{\}^+ { } +
3 LAR与LASSO
LAR虽然是求解LASSO的算法，但它得到的解的路径，在出现了某个系数要穿过0 0 0
我们在第1节算法的第4步中加入一个规则：
若一个非零系数又变为了0 0 0
在修正后，LAR就可以解任意LASSO问题，包括p ≫ N p\gg N p ≫ N
为什么会出现与LASSO解不同的情况？我们注意到，对于LASSO的active set B \mathcal{B} B
x ⋅ j ′ ( y − X β ^ ) = λ sign ( β ^ j ) x_{\cdot j}'(y-X\hat\beta) = \lambda \text{sign}(\hat\beta_j) x ⋅ j ′ ( y − X β ^ ) = λ sign ( β ^ j ) 而对于LAR的active set A \mathcal{A} A
x ⋅ j ′ ( y − X β ^ ) = γ s j x_{\cdot j}'(y-X\hat\beta) = \gamma s_j x ⋅ j ′ ( y − X β ^ ) = γ s j 其中s j s_j s j
在正常情况下，上面二者的右侧是相等的，也因此LAR的解就是LASSO的解。但是，当一个非零系数要穿过0 0 0 B \mathcal{B} B s j s_j s j 0 0 0
x ⋅ j ′ ( y − X β ^ ) ≤ λ x_{\cdot j}'(y-X\hat\beta) \leq \lambda x ⋅ j ′ ( y − X β ^ ) ≤ λ 这恰恰与A c \mathcal{A}^c A c A \mathcal{A} A
参考文献
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.
