第3章 - Least Angle Regression

我会慢慢补齐这一份笔记的。

首先定义一下数据矩阵($X_{(N\times p)}$)，每行一个数据点，响应向量（response）($\mathbf{y}_{(N\times 1)}$)，数据已经过标准化处理，即满足：

($\sum_{i=1}^N y_i = 0, \quad \sum_{i=1}^N x_{ij} = 0,  \quad \sum_{i=1}^N x_{ij}^2 = 1$)

并假设($X$)列向量线性独立。

一、LARS的计算步骤

LARS每次选择一个变量（predictor），加入到集合($\mathcal{A} \subset \{1,2,\cdots,p\}$)中。假设在第k次，已经选出了一个active变量集合($\mathcal{A}$)，同时已经有了一个对($\mathbf{y}$)的估计($\hat{\mathbf{y}}_{\mathcal{A}}$)。现在第k+1步修正估计($\hat{\mathbf{y}}$)：

($\begin{equation}\label{lars-update}
\hat{\mathbf{y}}_{\mathcal{A}+} = \hat{\mathbf{y}}_\mathcal{A} + \hat{\gamma} \mathbf{u}_\mathcal{A}
\end{equation}$)

其中($\mathbf{u}_{\mathcal{A}}$)和($\hat{\gamma}$)按如下计算：

(1) 定义数据和当前残差的correlation向量

($\hat{\mathbf{c}} = X^T(\mathbf{y} - \hat{\mathbf{y}}_{\mathcal{A}})$)

LARS的策略会使得当前active集合($\mathcal{A}$)中所有的predictor和残差的corr绝对值都相等，并且是最大的。也就是说

($\hat{C} = \max \limits_j \{| \hat{c}_j |\} \qquad \mathcal{A} = \{ j : |\hat{c}_j|=\hat{C} \}$)

（为什么是这样，后面会解释）

(2) 定义一个矩阵

($X_{\mathcal{A}} = \left( \cdots s_j \mathbf{x}_j \cdots \right)_{j\in \mathcal{A}}$)

其中($s_j = sign \{\hat{c}_j\} \in \{+1, -1\}$)

(3) 定义两个中间量

($\mathcal{G}_\mathcal{A} = X_{\mathcal{A}}^T X_{\mathcal{A}} \qquad A_\mathcal{A} = (\mathbf{1}_\mathcal{A}^T \mathcal{G}_\mathcal{A}^{-1} \mathbf{1}_\mathcal{A})^{-1/2}$)

其中($\mathbf{1}_\mathcal{A}$)是长度为($|\mathcal{A}|$)的全1向量。

(4) 得到搜索方向($\mathbf{u}_\mathcal{A}$)

($\mathbf{u}_\mathcal{A} = X_{\mathcal{A}} w_\mathcal{A}
\quad where \quad
w_\mathcal{A} = A_\mathcal{A} \mathcal{G}_\mathcal{A}^{-1} \mathbf{1}_\mathcal{A}$)

容易验证，($\mathbf{u}_\mathcal{A}$)是($\mathcal{A}$)中所有predictor的角平分线（夹角小于90度），并且是一个单位方向，即：

($X_\mathcal{A}^T \mathbf{u}_\mathcal{A} = A_\mathcal{A} \mathbf{1}_\mathcal{A}
\quad and \quad \|\mathbf{u}_\mathcal{A}\|^2 = 1$)

(5) 得到搜索步长($\hat{\gamma}$)

($\hat{\gamma} = \min^+_{j\in \mathcal{A}^c} \left\{ \frac{\hat{C}-\hat{c}_j}{A_\mathcal{A}-a_j}, \frac{\hat{C}+\hat{c}_j}{A_\mathcal{A}+a_j} \right\}$)

其中($a_j = \mathbf{x}_j^T\mathbf{u}_\mathcal{A}$)。($\min^+_{j \in \mathcal{A}^c}$)表示在等号右侧的两个值中选正的那个，然后在inactive集合($\mathcal{A}^c$)中选择一个使得这个值最小的predictor，加入active集合($\mathcal{A}_+$)中。

至此，LARS的基本算法就描述完整了。另外第一步的选择其实没有什么特殊，令第0步的估计($\hat{\mathbf{y}}_{\mathcal{A}_0} = \vec{0}$)，其中($\mathcal{A}_0$)里有一个和现有残差|corr|最大的($\mathbf{x}_{j1}$)，然后计算可知，下一步的搜索方向就是这个($\mathbf{x}_{j1}$)的方向。

二、LARS的性质

下面我们看一下选择的角平分线方向($\mathbf{u}_\mathcal{A}$)和步长($\hat{\gamma}$)有怎么样的性质。定义

($\mathbf{y}_{\mathcal{A}_+}(\gamma) \equiv \hat{\mathbf{y}}_\mathcal{A} + \gamma \mathbf{u}_\mathcal{A}$)

其中（根据定义）($\gamma > 0$)，于是对于任意一个predictor，它和 改进后的残差 的corr是这样的：

($c_j(\gamma) = \mathbf{x}_j^T (\mathbf{y} - \mathbf{y}_{\mathcal{A}_+}(\gamma))
= \hat{c}_j - \gamma \mathbf{x}_j^T \mathbf{u}_\mathcal{A} \equiv \hat{c}_j -\gamma a_j$)

1. 对于($j \in \mathcal{A}$)，我们发现：
若($\hat{c}_j > 0$)，则($s_j = +1$)，则($\mathbf{x}_j^T \mathbf{u}_\mathcal{A} = s_j \mathbf{x}_j^T \mathbf{u}_\mathcal{A} = A_\mathcal{A}$)，即($\hat{c}_j - \gamma \mathbf{x}_j^T \mathbf{u}_\mathcal{A} = \hat{C} - \gamma A_\mathcal{A}$)
若($\hat{c}_j < 0$)，则($s_j = -1$)，则($\mathbf{x}_j^T \mathbf{u}_\mathcal{A} = -s_j \mathbf{x}_j^T \mathbf{u}_\mathcal{A} = -A_\mathcal{A}$)，即($\hat{c}_j - \gamma \mathbf{x}_j^T \mathbf{u}_\mathcal{A} = -\hat{C} + \gamma A_\mathcal{A}$)

令($\hat{C} - \gamma A_\mathcal{A} > 0$)。综合得到，($|c_j(\gamma)| = \hat{C} - \gamma A_\mathcal{A}$)

也就是说，所有($j \in \mathcal{A}$)同新的残差方向的corr（的绝对值）仍然相等（减小量相等）。

2. 对于($j \in \mathcal{A}^c$)，我们令($|c_j(\gamma)| = \hat{C} - \gamma A_\mathcal{A}$)，就解出了上文第(5)步等式右边的表达式，我们选择最小的那个步长，并把这个($j$)加入新的($\mathcal{A}_+$)

至此我们看到，角平分线方向($\mathbf{u}_\mathcal{A}$)和步长($\hat{\gamma}$)的选择令现有active集合中所有的predictor和新残差的($|\operatorname{corr}_\mathcal{A}^{new}|$)保持一致，并且这个($|\operatorname{corr}_\mathcal{A}^{new}|$)也等于新加入的predictor和新残差的($|\operatorname{corr}_{\mathcal{A}_+}^{new}|$)。

读ESL的朋友会发现，上文对LAR算法的描述和书中的描述并不一致，书里面是这么写的：

1. Standardize the predictors to have mean zero and unit norm. Start with the residual ($r = y − \bar{y}$)，($ \beta_1 , \beta_2 , \dots , \beta_p = 0$).
2. Find the predictor ($x_j$) most correlated with ($r$).
3. Move ($\beta_j$) from 0 towards its least-squares coefficient ($x_j$) , ($r$) , until some other competitor ($x_k$) has as much correlation with the current residual as does ($x_j$) .
4. Move ($\beta_j$) and ($\beta_k$) in the direction defined by their joint least squares coefficient of the current residual on ($(x_j , x_k)$), until some other competitor ($x_l$) has as much correlation with the current residual.
5. Continue in this way until all ($p$) predictors have been entered. After ($\min(N − 1, p)$) steps, we arrive at the full least-squares solution.

第5步所谓『continue』指的是，接下来把($(\beta_j，\beta_k，\beta_l)$)向($(x_j，x_k，x_l)$)（以新残差为回归目标）的最小二乘方向移动。以此类推，直到所有的变量都被选入。

下面我们来看一下，为什么这两个算法描述是等价的。

首先定义一下符号。令($\tilde{X}_\mathcal{A}$)表示所有已经选入的predictor组成的矩阵，即

($\tilde{X}_\mathcal{A} = X_\mathcal{A} S_\mathcal{A},
\quad
S_\mathcal{A} = \left(
  \begin{array}{cccc}
    s_1  \\
    & s_2  \\
    && \ddots \\
    &&& s_{|\mathcal{A}|}
  \end{array}
\right)
$)

定义($\Delta\mathbf{y}_\mathcal{A} \equiv \mathbf{y} - \hat{\mathbf{y}}_{\mathcal{A}}$)，我们知道其最小二乘的解：

($\begin{align}
minimize \frac{1}{2} \| \Delta\mathbf{y}_\mathcal{A} - 
\tilde{X}_\mathcal{A} \pmb{\beta}_\mathcal{A} \|_2^2 
\quad &\Rightarrow \quad 
\hat{\pmb{\beta}}_\mathcal{A} = 
(\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1} \tilde{X}_\mathcal{A}^T \Delta\mathbf{y}_{\mathcal{A}} 
= (\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1} \hat{\mathbf{c}}_\mathcal{A} \\
&\Rightarrow \quad
\tilde{X}_\mathcal{A}\hat{\pmb{\beta}}_\mathcal{A}
= \hat{C}\tilde{X}_\mathcal{A}(\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1} \mathbf{s}_\mathcal{A}
\end{align}$)

其中($\mathbf{s}_\mathcal{A} = [\cdots s_j \cdots ]_{j\in\mathcal{A}}^T$)

接下来我们对($\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A}$)做SVD分解：

($
\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A} 
= U\Sigma V^T \quad \Rightarrow \quad 
(\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1}
= V\Sigma^{-1}U^T
$)

两边乘($S_\mathcal{A}$)，

($
X_\mathcal{A}^T X_\mathcal{A} =
S_\mathcal{A}\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A}S_\mathcal{A}
= S_\mathcal{A}U\Sigma V^T S_\mathcal{A}
$)

可见($S_\mathcal{A}U$)和($S_\mathcal{A}V$)是正交矩阵，所以上式是($X_\mathcal{A}^T X_\mathcal{A}$)的SVD分解，因此，

($\mathcal{G}_\mathcal{A}^{-1} = (X_\mathcal{A}^T X_\mathcal{A})^{-1}
= S_\mathcal{A}V\Sigma^{-1} U^T S_\mathcal{A}
= S_\mathcal{A} (\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1} S_\mathcal{A}
$)

由此得到，

($\mathbf{u}_\mathcal{A} = X_{\mathcal{A}} w_\mathcal{A}
=
\tilde{X}_{\mathcal{A}} S_\mathcal{A} A_\mathcal{A} \mathcal{G}_\mathcal{A}^{-1} \mathbf{1}_\mathcal{A}
= A_\mathcal{A} \tilde{X}_\mathcal{A}(\tilde{X}_\mathcal{A}^T\tilde{X}_\mathcal{A})^{-1} \mathbf{s}_\mathcal{A}
$)

现在就比较清楚了，LAR系数更新的方向确实就是（以当前残差为目标）最小二乘的方向。而这个方向上的步长，如本节2.的描述，就是Hastie所言『until some other competitor ($x_l$) has as much correlation with the current residual』。

三、LARS的Lasso修正

首先上一张图，来自ESL figure 3.15：
LAR和Lasso求解路径对比

左边是LAR各变量的求解路径，($\pmb{\beta}$)各维依次被加入active集合（从0变化到非0）。右边是Lasso的。可见两张图几乎一样，区别起始于大约横坐标=18的地方，深蓝色的那条线重新变成了0。下面我们来仔细分析一下这里究竟发生了什么。

首先来看Lasso的loss function：

($R(\pmb{\beta}) = \frac{1}{2} \| \mathbf{y} - \mathbf{X}\pmb{\beta}\|_2^2 + \lambda \| \pmb{\beta} \|_1$)

虽然L1-norm是不可导的，但对于已经加入active set ($\mathcal{A}$)的($\beta_j \neq 0$)来说，我们可以放心地求导：

($\frac{\partial R}{\partial \pmb{\beta}} = -\mathbf{X}^T\mathbf{y} + \mathbf{X}^T\mathbf{X}\pmb{\beta} + \lambda \operatorname{sign} (\pmb{\beta}), \quad \forall j \in \mathcal{A}$)

令上式=0，解出($\pmb{\beta}$)的估计($\hat{\pmb{\beta}}$)。于是我们有，

($\hat{c}_j \equiv \mathbf{x}_j^T (\mathbf{y} - \hat{\mathbf{y}}) = \lambda \operatorname{sign}(\hat{\beta_j}), \quad \forall j \in \mathcal{A}$)

其中($\hat{\mathbf{y}}=\mathbf{X}\pmb{\beta}$)。由于($\lambda>0$)，我们有

($\begin{equation}\label{lasso-sign-equation}
\operatorname{sign}(\hat{\beta}_j) = \operatorname{sign}(\hat{c}_j) = s_j
\end{equation}$)

我们记住Lasso的这一结论，后面会用到。

接下来我们来看LARS的求解。假设我们现在已经有了一个active集合($\mathcal{A}$)，一个和Lasso估计($\hat{\mathbf{y}} = \mathbf{X}\hat{\pmb{\beta}}$)对应的LARS估计($\hat{\mathbf{y}}_\mathcal{A} = \mathbf{X}\hat{\pmb{\beta}}_\mathcal{A}$)，下面做一步LARS更新：

($\begin{align}
\hat{\mathbf{y}}_{\mathcal{A}+}(\gamma) 
&= \hat{\mathbf{y}}_\mathcal{A} + \gamma\mathbf{u}_\mathcal{A}  \\
&= \mathbf{X}\hat{\pmb{\beta}}_\mathcal{A} + \gamma \mathbf{X}_\mathcal{A} \mathbf{w}_\mathcal{A} \\
&= \mathbf{X} (\hat{\pmb{\beta}}_\mathcal{A} + \gamma \hat{\mathbf{d}} )
\end{align}$)

这里的($\hat{\mathbf{d}}$)是新定义的一个向量：

($
\begin{equation}
\hat{d}_j=\begin{cases}
s_j w_{\mathcal{A} j} \quad  \text{for}  \quad j \in \mathcal{A}\\
0  \quad \text{otherwise}
\end{cases}
\end{equation}
$)

这么一来，我们得到新的coefficients：

($\beta_{\mathcal{A}_+j}(\gamma) = \hat{\beta}_{\mathcal{A}j} + \gamma \hat{d}_j$)

注意到从($\gamma_j = -\hat{\beta}_{\mathcal{A}j} / \hat{d}_j$)开始，($\beta_{\mathcal{A}_+j}(\gamma)$)的符号就改变了。由于LARS每次总是取满足条件的最小的步长，我们定义：

($
\begin{equation}
\widetilde{\gamma}=\begin{cases}
\min\{\gamma_j\} \quad  存在\gamma_j > 0\\
+\infty  \quad \text{otherwise}
\end{cases}
\end{equation}
$)

若是LARS更新步长($\hat{\gamma} > \widetilde{\gamma}$)，那么这次更新过后，($\beta_j$)的符号改变了，而($c_j$)的符号不变（($c_j = \hat{C}-\gamma A_\mathcal{A}$) or ($c_j = -\hat{C}+\gamma A_\mathcal{A}$)），如此就和Lasso的要求（两者的符号始终一致）矛盾。这就是LARS和Lasso求解路径的差异所在（图中($\beta$)重新变为0就是上文中所说的($\beta$)符号改变的情况）。

解决方法很简单，就是对($\hat{\gamma} > \widetilde{\gamma}$)的情况特殊处理，并把这个predictor从active集合中挪走，即：

($\hat{\mathbf{y}}_{\mathcal{A}+} = \hat{\mathbf{y}}_\mathcal{A} + \widetilde{\gamma} \mathbf{u}_\mathcal{A} \quad and \quad \mathcal{A}_+ = \mathcal{A} - \{\widetilde{j}\}$)

四、LARS和Lasso的关系

那么LARS和Lasso究竟是怎么个关系呢，它们俩的求解路径为什么这么像？LARS的搜索方向($\mathbf{u}_\mathcal{A} = X_{\mathcal{A}} w_\mathcal{A}$)和Lasso有什么关联？现在我们就来仔细考察一下。

先来回顾一下Lasso的优化目标：

($minimize \| \mathbf{y} - \mathbf{X}\pmb{\beta} \|_2^2 $)
($st \quad \|\pmb{\beta}\|_1 \leq t$)

我们希望通过KKT条件来分析最优解的性质，但是这个L1-regularization真的很碍眼：它不能求导。于是我们把($\pmb{\beta}$)拆成两部分：

($\pmb{\beta} = \pmb{\beta}_+ - \pmb{\beta}_- \qquad 
\pmb{\beta}_+ \succeq \mathbf{0}, \quad \pmb{\beta}_- \succeq \mathbf{0}$)

其中($\succeq$)表示向量逐元素大于等于，即在($\mathbf{R}_+$)这个正常锥定义下的广义不等式。

然后我们就可以重写优化目标如下：

($minimize \quad \frac{1}{2} \sum_{i=1}^N \left [ y_i - (\sum_{j=1}^p x_{ij} \beta_{j+} - \sum_1^p x_{ij} \beta_{j-} ) \right ]^2$)
($st \quad
-\beta_{j+} \leq 0, \quad
-\beta_{j-} \leq 0, \quad
\sum_{j=1}^p (\beta_{j+} + \beta_{j-}) \leq t$)

为了让大家看得清楚明白不留死角，我把完整的计算过程写出来：

($\begin{align}
\frac{\partial \left\{ \frac{1}{2} \sum_{i=1}^N \big [ y_i - (\sum_{j=1}^p x_{ij} \beta_{j+} - \sum_{j=1}^p x_{ij} \beta_{j-} ) \big ]^2 \right \}}{\partial \beta_{j+}} &= \sum_{i=1}^N \left \{ \big [ y_i - (\sum_{j=1}^p x_{ij} \beta_{j+} - \sum_{j=1}^p x_{ij} \beta_{j-} ) \big ] \cdot (-x_{ij}) \right \} \\
\frac{\partial \left\{ \frac{1}{2} \sum_{i=1}^N \big [ y_i - (\sum_{j=1}^p x_{ij} \beta_{j+} - \sum_{j=1}^p x_{ij} \beta_{j-} ) \big ]^2 \right \}}{\partial \beta_{j-}} &= \sum_{i=1}^N \left\{ \big [ y_i - (\sum_{j=1}^p x_{ij} \beta_{j+} - \sum_{j=1}^p x_{ij} \beta_{j-} ) \big ] \cdot x_{ij} \right\} \\
\frac{\partial \left\{ -\beta_{j+} \right\}}{\partial \beta_{j+}} &= -1 \\
\frac{\partial \left\{ -\beta_{j-} \right\}}{\partial \beta_{j-}} &= -1 \\
\frac{\partial \left\{ \sum_{j=1}^p (\beta_{j+} + \beta_{j-}) -t \right\}}{\partial \beta_{j+}} &= 1 \\
\frac{\partial \left\{ \sum_{j=1}^p (\beta_{j+} + \beta_{j-}) -t \right\}}{\partial \beta_{j-}} &= 1
\end{align}$)

对可行性条件分别引入Lagrange乘子($\lambda_{j+} \ge 0$)、($\lambda_{j-}\ge 0$)、($\lambda\ge 0$)，我们得到Lasso的KKT最优条件如下：

($\begin{align}
-\mathbf{x}_j^T\mathbf{r} + \lambda - \lambda_{j+} &= 0 \\
\mathbf{x}_j^T\mathbf{r} + \lambda - \lambda_{j-} &= 0 \\
\lambda_{j+}\beta_{j+} &= 0 \\
\lambda_{j-}\beta_{j-} &= 0 \\
\lambda \left[ \sum_{j=1}^p (\beta_{j+} + \beta_{j-}) -t \right] &= 0
\end{align}$)

其中($\mathbf{r} = \mathbf{y} - \left( \sum_{j=1}^p \mathbf{x}_j \beta_{j+} -  \sum_{j=1}^p \mathbf{x}_j \beta_{j-} \right)$)

由以上等式，我们可以作如下推导：

($\begin{align}
\lambda=0 
&\Rightarrow \lambda_{j-} = \mathbf{x}_j^T\mathbf{r} = -\lambda_{j+} \\
&\Rightarrow \mathbf{x}_j^T\mathbf{r} = 0
\end{align}$)

这种情况就退化为不带regularization的一般OLS（回忆least square求导等于0就是上式）。于是我们就只要考虑($\lambda >0$)的情况：

($\begin{align}
\beta_{j+} > 0, \lambda > 0 
&\Rightarrow \lambda_{j+} = 0 \\
&\Rightarrow \mathbf{x}_j^T\mathbf{r} = \lambda > 0 \\
&\Rightarrow \lambda_{j-} = \mathbf{x}_j^T\mathbf{r} + \lambda > 0 \\
&\Rightarrow \beta_{j-} = 0
\end{align}$)

同理可得：

($\begin{align}
\beta_{j-} > 0, \lambda > 0 
&\Rightarrow \mathbf{x}_j^T\mathbf{r} = -\lambda < 0\\
&\Rightarrow \lambda_{j+} = 0
\end{align}$)

也就是说，($\mathbf{x}_j$)和残差的correlation总是同其对应的($\beta_j$)符号一致。如果大家看Efron的原始论文，(5.19)式下面有一句话说『... implies we can drop the second condition in (5.19) ...』，完整的证明就如上文所述。Efron段位太高，他看一眼就知道的事情，本人愚钝，参考了杨灿老师的pdf才勉强证出。

接下来我们试图把LAR和Lasso归入一个框架体系中。LAR中的active集合($\mathcal{A}$)在这里就是所有不为0的($\beta_j$)：

($\mathcal{A} = \{ j: \beta_{j+} > 0 || \beta_{j-} > 0\}$)

令($\mathbf{S}_\mathcal{A} = diag(s_1, s_2, \dots, s_{|\mathcal{A}|})$)，即对角线元素是($sign(\beta_j)$)的对角矩阵。根据之前的推导，我们可以得到：

($\mathbf{S}_\mathcal{A} \mathbf{X}_\mathcal{A}^T \big[ \mathbf{y} - \mathbf{S}_\mathcal{A} \mathbf{X}_\mathcal{A} \hat{\pmb{\beta}}_\mathcal{A}(t) \big] = \lambda \mathbf{S}_\mathcal{A} \vec{1}$)

分别对两个不同的正则约束($t_1 < t_2$)（相应的Lagrange乘子($\lambda_1$)和($\lambda_2$)），并在上式两边左乘($\mathbf{S}_\mathcal{A}$)，我们得到：

($\begin{align}
\mathbf{X}_\mathcal{A}^T\mathbf{X}_\mathcal{A}\mathbf{S}_\mathcal{A} ( \hat{\pmb{\beta}}_\mathcal{A}(t_2) - \hat{\pmb{\beta}}_\mathcal{A}(t_1)) &= (\lambda_1 - \lambda_2) \vec{1} \\
\Rightarrow \hat{\pmb{\beta}}_\mathcal{A}(t_2) - \hat{\pmb{\beta}}_\mathcal{A}(t_1) &= (\lambda_1 - \lambda_2) \mathbf{S}_\mathcal{A} \mathbf{\mathcal{G}}_\mathcal{A}^{-1} \vec{1}
\end{align}$)

然后我们从L1-regularization的定义中引入
($\| \hat{\pmb{\beta}} \|_1 = \sum_\mathcal{A} s_j \hat{\beta}_j = t$)
($\Rightarrow \mathbf{s}_\mathcal{A}^T[ \hat{\pmb{\beta}}_\mathcal{A}(t_2) -\hat{\pmb{\beta}}_\mathcal{A}(t_1)] = t_2 - t_1$)

其中($\mathbf{s}_\mathcal{A}$)是一个列向量，对应位置元素为($s_j$)。注意到我这里偷偷把正则约束的($\leq$)换成了($=$)，这是因为，只要($t < \| \hat{\pmb{\beta}}_{ols} \|_1$)，也就是当约束起作用的时候，将不等号替换为等号得到的约束效果是一致的。

于是综合得到：

($t_2 - t_1 = (\lambda_1 - \lambda_2)\mathbf{s}_\mathcal{A}^T\mathbf{S}_\mathcal{A}\mathbf{\mathcal{G}}_\mathcal{A}^{-1} \vec{1} = (\lambda_1 - \lambda_2) \vec{1}^T \mathbf{\mathcal{G}}_\mathcal{A} \vec{1} = (\lambda_1 -\lambda_2) A_{\mathcal{A}}^{-2}$)

($\hat{\pmb{\beta}}_\mathcal{A}(t_2) - \hat{\pmb{\beta}}_\mathcal{A}(t_1) 
= \mathbf{S}_\mathcal{A} A_\mathcal{A}^2 (t_2 - t_1) \mathbf{\mathcal{G}}_\mathcal{A}^{-1} \vec{1} = \mathbf{S}_\mathcal{A} A_\mathcal{A} (t_2-t_1) w_\mathcal{A}$)

对比之前的LARS的求解过程，两者更新($\pmb{\beta}$)的方向是一致的。因此，LARS可以认为是Lasso不断放松正则约束($t$)的过程。严格的证明还需要一些涉及细节的推论的支持，这部分就请拜读Efron教授的原文了。

五、(forward) stepwise & stagewise regression

下面来追溯一下LARS的源头。

我们找线性回归的稀疏解，并非一开始就有那么高级的LARS。最初的想法很简单，比如说一种叫做forward stepwise的方法，依照如下步骤：

1. 从inactive集合中找到一个变量($\mathbf{x}_{\hat{j}}$)：($\hat{j} = \operatorname{argmax}_{j\in\mathcal{A^c}} \left| \hat{c}_j \right|$)，其中($ \hat{c}_j \equiv \mathbf{x}_j^T (\mathbf{y}-\hat{\mathbf{y}}_{\mathcal{A}})$)

2. 定义($\Delta\mathbf{y}_\mathcal{A} \equiv \mathbf{y} - \hat{\mathbf{y}}_{\mathcal{A}}$)，用最小二乘计算其系数：

($minimize \frac{1}{2} \| \Delta\mathbf{y}_\mathcal{A} - \beta_{\hat{j}} \mathbf{x}_{\hat{j}}\|_2^2 \quad \Rightarrow \quad 
\hat{\beta}_{\hat{j}} = \Delta\mathbf{y}_{\mathcal{A}}^T \mathbf{x}_{\hat{j}} = \hat{c}_\hat{j}$)

（根据本文开头的标准化处理，($\mathbf{x}_j^T \mathbf{x}_j = 1$)）

3. 把($\mathbf{x}_{\hat{j}}$)加入active集合($\mathcal{A}$)，然后重复第1步。

我们注意到第2步计算得到的($\hat{\beta}_{\hat{j}}$)使得新加入的变量和新的残差($\Delta\mathbf{y}_{\mathcal{A}^+} \equiv \Delta\mathbf{y}_\mathcal{A} - \hat{\beta}_{\hat{j}} \mathbf{x}_{\hat{j}}$)正交了（可以自己带入算一下）：

($\mathbf{x}_{\hat{j}}^T \Delta\mathbf{y}_{\mathcal{A}_+}= 0$)

这意味着什么呢？由于下一个变量($\mathbf{x}_{\hat{j}_+}$)是向($\Delta\mathbf{y}_{\mathcal{A}_+}$)回归，那么($\mathbf{x}_{\hat{j}_+}$)和($\mathbf{x}_{\hat{j}}$)就会比较接近正交。Efron教授指出，stepwise的步长太激进了，从而可能导致一些和($\mathbf{x}_{\hat{j}}$)有较大corr的变量没有机会被选中。改进方法很简单：我们在每个步长前面乘上一个很小的系数：

($\hat{\mathbf{y}} \gets \hat{\mathbf{y}} + \epsilon \operatorname{sign}(\hat{c}_\hat{j}) \mathbf{x}_\hat{j}$)

Efron教授把这个东西叫做forward stagewise regression。为什么叫『stage』呢？看下面这张图（取自Efron B, Hastie T, Johnstone I, et al. Least angle regression[J]. The Annals of statistics, 2004, 32(2): 407-499.）：
LAR和stagewise regression

如果当前的估计为($\hat{\mathbf{\mu}}_1$)，则当前残差($\bar{\mathbf{y}}_2 - \hat{\mathbf{\mu}}_1$)和($\mathbf{x}_1$)、($\mathbf{x}_2$)的corr相等，我们随便挑选一个，比如($\mathbf{x}_1$)，作为搜索方向，移动一小步。这时候新的残差就和($\mathbf{x}_2$)有了更大的corr，然后我们选择($\mathbf{x}_2$)作为搜索方向，移动一小步。注意这里的图解有些误导，竖直的step其实应该和($\mathbf{x}_2$)的方向平行。如此反复，就跟阶梯（stage）一样。

我们后面将看到严格的证明：当($\epsilon \rightarrow 0$)的时候，stagewise和（经改进的）LARS是一致的。

六、LARS的stagewise修正

Efron本来觉得，当($\epsilon \rightarrow 0$)的时候，stagewise就是LARS。从上一节的图像解释看，似乎的确如此。事实上，未做修正的LARS只在($p=2$)（也就是数据是2维）的情况下和stagewise一致，对于高维情况，需要做一个简单的修正。下面我们就来看看造成差异的原因和修正方法。

假设对某次stagewise而言，已经有了一个估计($\hat{\mathbf{y}}$)，然后在这个估计的基础上移动了N个step，每个step都乘上了很小的($\epsilon$)，即($\epsilon s_j \mathbf{x}_j$)。我们把($N_j$)定义为这期间某个($\mathbf{x}_j$)被选为step的次数。容易证明，当($\epsilon$)足够小的情况下，在这N步期间inactive集合中的元素不可能被选出，即($N_j = 0, \forall j \in \mathcal{A}^c$)（本节末尾会给出证明）

下面我们定义($\pmb{\rho} \equiv (N_1, N_2, \dots , N_m) / N$)，并令($\pmb{\rho}_\mathcal{A}$)表示active集合构成的($\pmb{\rho}$)向量，得到新的估计：

($\hat{\mathbf{y}}_+ = \hat{\mathbf{y}} + N\epsilon \mathbf{X}_{\mathcal{A}} \pmb{\rho}_\mathcal{A}$)

对比之前得出的LARS的更新策略：

($\hat{\mathbf{y}}_{\mathcal{A}+} 
= \hat{\mathbf{y}}_\mathcal{A} + \hat{\gamma} \mathbf{u}_\mathcal{A}
= \hat{\mathbf{y}}_\mathcal{A} + \hat{\gamma} \mathbf{X}_\mathcal{A} w_\mathcal{A} \quad where \quad 
w_\mathcal{A} = A_\mathcal{A} \mathcal{G}_\mathcal{A}^{-1} \vec{1}_\mathcal{A}$)

我们发现，由于($\pmb{\rho}_\mathcal{A} \succeq \vec{0}$)，当($w_\mathcal{A} <0$)时，LARS和stagewise的搜索路径就不一致了。解决方法就是把($\mathbf{u}_\mathcal{A}$)扔到如下这个凸锥里去：

($\mathcal{C}_\mathcal{A} = \left \{ \mathbf{v} = \sum\limits_{j\in \mathcal{A}} s_j \mathbf{x}_j \rho_j, \quad \rho_j \geq 0 \right \}$)

于是LARS的stagewise修正如下：若($\mathbf{u}_\mathcal{A} \in \mathcal{C}_\mathcal{A}$)，则执行正常的LARS更新；反之，找到($\mathbf{u}_\mathcal{A}$)在($\mathcal{C}_\mathcal{A}$)里的投影，即这个凸锥中离($\mathbf{u}_\mathcal{A}$)最近的点($\mathbf{u}_\mathcal{A}^{proj}$)，在这个方向上更新。

========================割：证明本节开头那个结论====================

下面我们证明：

对于足够小的($\epsilon$)，在N步之内，($N_j = 0, \forall j \in \mathcal{A}^c$)

定义($\gamma \equiv N\epsilon$)，前一次估计为($\hat{\mathbf{y}}$)，N步之后的估计为($\hat{\mathbf{y}}(N)$)，定义($\hat{\mathbf{c}}(N) \equiv \mathbf{X}^T(\mathbf{y} - \hat{\mathbf{y}}(N))$)，有：

($\|\hat{\mathbf{y}}(N) - \hat{\mathbf{y}}\| \leq \gamma
\Rightarrow |\hat{c}_j(N) - \hat{c}_j| = |\mathbf{x}_j^T(\hat{\mathbf{y}}(N) - \hat{\mathbf{y}})|
\leq \|\mathbf{x}_j\| \cdot \| \hat{\mathbf{y}}(N) - \hat{\mathbf{y}} \| \leq \gamma
$)

可见，令($\Delta_{\mathcal{A}^c} = \hat{C} - \max_{\mathcal{A}^c}\{ \hat{c}_j \}$)，当($\gamma < \frac{1}{2} \Delta_{\mathcal{A}^c}$)时，inactive集合($\mathcal{A}^c$)里的元素总是无法被选出。（这是由于N次更新前后($c_j$)的变化量总是小于一半($\Delta_{\mathcal{A}^c}$)，那么极端情况，inactive集合中拥有最大corr的($\mathbf{x}_{\mathcal{A}^c}$)增加了一半($\Delta_{\mathcal{A}^c}$)，active集合中的($\mathbf{x}_{\mathcal{A}}$)减少了一半($\Delta_{\mathcal{A}^c}$)，前者还是干不过后者，无法被选出。）

七、stagewise和LAR的关系

对于stagewise regression，其实我们还有一件事没有讲，就是在实际应用中，某一个predictor更新次数($N_j$)的问题。既然我们的目的是要减小残差，那么active集合中的predictor自然要尽可能地朝残差减小的方向移动，并满足一些限制。也就是求解如下问题（($\rho_j = N_j / N$)）：

($minimize \quad \frac{1}{2} \| \mathbf{r} - N \epsilon \mathbf{X}_\mathcal{A} 
\pmb{\rho}_\mathcal{A}  \|_2^2
\quad
st \quad
\pmb{\rho} \succeq \vec{0}, \quad
\pmb{\rho}^T \vec{1} = 1$)

照例来考察KKT最优条件：

($\begin{align}
\frac{\partial \left\{ \frac{1}{2} \| \mathbf{r} - N \epsilon X_\mathcal{A}\pmb{\rho}_\mathcal{A} \|_2^2 \right\}}{\partial \pmb{\rho}_\mathcal{A}} 
&= -N\epsilon X_\mathcal{A}^T \mathbf{r} + N^2\epsilon^2 X_\mathcal{A}^T X_\mathcal{A} \pmb{\rho}_\mathcal{A} \\
\frac{\partial \left\{ -\pmb{\rho}_\mathcal{A} \right \}}{\partial \pmb{\rho}_\mathcal{A}} &= -I \\
\frac{\partial \left\{ \pmb{\rho}_\mathcal{A}^T \vec{1} -1 \right\}}{\partial \pmb{\rho}_\mathcal{A}} &= \vec{1}
\end{align}$)

分别引入Lagrange乘子($\pmb{\gamma} \succeq \vec{0}$), ($\lambda \ge 0$)，得到KKT条件为：

($\begin{align}
 -N\epsilon X_\mathcal{A}^T \mathbf{r} + N^2\epsilon^2 X_\mathcal{A}^T X_\mathcal{A} \pmb{\rho}_\mathcal{A} - \pmb{\gamma} + \lambda \vec{1}
&= \vec{0} \\
\pmb{\rho}_\mathcal{A}^T\vec{1} - 1 &= 0 \\
\lambda &\ge 0 \\
\pmb{\rho}_\mathcal{A} &\succeq \vec{0} \\
\pmb{\gamma} &\succeq \vec{0} \\
\pmb{\gamma}^T \pmb{\rho}_\mathcal{A} &= 0 \\
\end{align}$)
注意到联立最后三条可知，($\gamma_j \rho_j = 0$)。

KKT的第一个条件告诉我们，($-N\epsilon X_\mathcal{A} (\mathbf{r} - N\epsilon X_\mathcal{A} \pmb{\rho} ) = \pmb{\gamma} + \lambda \vec{1}$)，也就是说，active集合中的predictor与残差的corr是相等的，这也是LAR算法的性质。

Hestie在（Hastie T, Taylor J, Tibshirani R, et al. Forward stagewise regression and the monotone lasso[J]. Electronic Journal of Statistics, 2007, 1: 1-29.）中证明了，上述优化问题和非负最小二乘：

($minimize \frac{1}{2} \| \mathbf{r} - X_\mathcal{A} \theta \|_2^ \quad st \quad \theta \succeq \vec{0}$)

给出相同的搜索方向，这也是stagewise regression和lasso一个变种『monotone lasso』相联系的重要一环。

至此，我们看到Stagewise / LAR / Lasso的搜索方向本质上是一致的，只是各自有细微不同的限制条件，我先抄一段杨灿老师ppt中的summary，回头再来仔细补：
LAR uses least squares directions in the active set of variables.
Lasso uses least square directions; if a variable crosses zero, it is removed from the active set.
Boosting uses non-negative least squares directions in the active set.