3.8 更多关于套索回归和类似的路径算法

最小角回归算法的发表（Efron et al., 2004）开启了一系列不同问题的正则化系数变化路径拟合的算法的研究。另外， $L_{1}$ 正则项的使用引导了信号处理中 压缩感知（compressed sense） 领域的发展。（Donoho, 2006a; Candes, 2006）本节从最小角回归算法的前身引入，并介绍一些相关的修改版本和其他的路径算法。

3.8.1 增量前向分段回归

增量前向分段回归（Incremental Forward Stagewise Regression） 是另一种类似于最小角回归算法的方法。有意思的是，作为新的线性模型算法的最小角回归，灵感来自于对灵活的非线性回归过程（提升方法）的研究。作者的同事 Brad Efron 在阅读了本书第十六章的前向分段算法 16.1 后¹意识到，在线性模型中可以直接构建如图 3.10 中的分段线性的套索系数路径。他也从此提出了第 3.4.4 节中的最小角回归算法，和本节的增量版本前向分段回归。

将第 16.1 节中的前向分段提升算法 16.1 应用在线性回归的场景中，不断地更新与当前残差相关性最高的变量的系数（改变幅度为一个小值 $ϵ$ ），可得到系数的变化曲线。算法 3.4 描述了计算细节。

算法 3.4：增量前向分段回归 $F S_{ϵ}$

将 $y$ 作为初始化残差 $r$ ，系数 $β_{1}, β_{2}, \dots, β_{p} = 0$ 。所有的输入变量均被标准化为均值 0 和范式（方差) 1。
找到与 $r$ 相关性最高的自变量 $x_{j}$ 。
更新对应的系数： $β_{j} \leftarrow β_{j} + δ_{j}$ 。其中 $δ_{j} = ϵ \cdot sign [⟨ x_{j}, r ⟩]$ ， $ϵ$ 为一个小的改变幅度。更新残差： $r \leftarrow r - δ_{j} x_{j}$
重复步骤 2 和 3，直到残差与所有的自变量都没有相关性。

**图 3.19**：前列腺癌症数据的系数曲线。左图为步幅为 $\epsilon=0.01$ 的增量前向分段回归。右图为通过 $\epsilon\Rightarrow 0$ 得出的无穷小量 $FS_0$。曲线通过对最小角回归算法 3.2 的修改版本 3.2b 计算得出。在此例中，$FS_0$ 的曲线是单调的，因此与套索回归和最小角回归的结果相同。 — **图 3.19**：前列腺癌症数据的系数曲线。左图为步幅为 $ϵ = 0.01$ 的增量前向分段回归。右图为通过 $ϵ \Rightarrow 0$ 得出的无穷小量 $F S_{0}$ 。曲线通过对最小角回归算法 3.2 的修改版本 3.2b 计算得出。在此例中， $F S_{0}$ 的曲线是单调的，因此与套索回归和最小角回归的结果相同。

图 3.19（左图）为在前列腺癌症数据中幅度为 $ϵ = 0.01$ 的算法结果。若算法中的 $δ_{j} = ⟨ x_{j}, r ⟩$ ，则与第 3.3.3 节中的前向分段回归过程完全一致。

我们主要感兴趣的是当改变幅度 $ϵ$ 很小时的性质。如果让 $ϵ \to 0$ ，则得到了图 3.19 的右图，其与图 3.10中的套索回归路径完全一致。我们将这个极限下的计算过程称为 无穷小量前向分段回归（infinitesimal forward stagewise regression），或 $F S_{0}$ 。这种过程在非线性自适应的方法中，例如第十和十六章中的提升方法，起着重要的作用，同时也是增量前向分段回归中最易于理论分析的版本。因其与提升方法的关联，Bühlmann and Hothorn (2007) 也将这个过程称为“L2 提升（L2boost）”。

Efron 原认为最小角回归的算法 3.2 可以实现 $F S_{0}$ 。使每个自变量保持与因变量一致的相关程度的同时，不断地逐渐更新相应的系数。然而他之后意识到最小角回归中对自变量的最小二乘拟合对某个系数的更新可能与对应自变量与因变量的相关方向（sign）相反，但这与算法 3.4 是不一致的。下面对最小角回归算法的修正可以使其完全实现 $F S_{0}$ ：

算法 3.2b：最小角回归的 $F S_{0}$ 修正

3.1b 计算含约束条件的最小回归估计，作为移动方向： $min_{b} ‖ r - X_{A} b ‖_{2}^{2} subject to b_{j} s_{j} \geq 0, j \in A$ 其中 $s_{j}$ 为 $⟨ x_{j}, r ⟩$ 的符号。

这个修正保证了系数的符号总是和相应输入变量与输出变量的相关性符号一致。可以证明（Hastie et al., 2007），在与输出变量相关性最大的一组输入变量上，这种方法是无穷小量更新幅度时的最优“更新步数（update turns）”。与套索回归类似， $F S_{0}$ 的完整系数路径可以通过（修正的）最小角回归算法快速的解出。

算法 3.2、3.2a、和 3.2b 的差别，导致当最小角回归算法的路径为如同图 3.19 中的单调递增或单调递减，则最小角回归、套索回归、和 $F S_{0}$ 会得到一致的路径。当路径不单调但不穿越 0 点，则最小角回归与套索回归一致。

$F S_{0}$ 是与最小角回归和套索回归相似又不同的方法，但却很难将其同样表达为对某个准则的最优化问题，因为 $F S_{0}$ 的系数路径实际上是一个微分方程的解。套索回归可理解为在固定的系数向量 $β$ 改变的 $L_{1}$ 范式下，寻找最小化残差平方和的变动方向； $F S_{0}$ 为在固定的系数变动路径的 $L_{1}$ 弧长下最小化，因而会阻止系数路径频繁地改变移动方向。

$F S_{0}$ 的约束性比套索回归更强，实际上可以看成是单调版本的套索回归。图 16.3 给出了一个比较极端的例子。 $F S_{0}$ 在 $p ≫ N$ 时可能仍有效，其系数曲线比套索回归更平滑，因此方差更小。更多 $F S_{0}$ 的内容可参考第 16.2.3 节以及 Hastie et al., 2007。图 3.16中的 $F S_{0}$ 系数曲线与套索回归非常相似。

3.8.2 分段线性路径算法

最小角回归算法利用了套索回归解的系数路径呈分段线性的特征，这也启发了其他正则问题中类似的路径算法。假设要求解：

\begin{matrix} (3.76) & \hat{β} (λ) = \underset{β}{\arg min} [R (β) + λ J (β)] \end{matrix}

\begin{matrix} (3.77) & R (β) = \sum_{i = 1}^{N} L (y_{i}, β_{0} + \sum_{j = 1}^{p} x_{i j} β_{j}) \end{matrix}

其中的损失函数 L 和惩罚项 J 均为凸函数。则解路径 $\hat{β} λ$ 为分段线性的充分条件为（Rosset and Zhu, 2007）：

R 为 $β$ 的二次或分段二次函数；
J 为 $β$ 的分段线性函数。

这也大体上意味着可以快速地计算解路径。例如，平方误差或绝对值误差损失函数，休伯化（Huberized）损失函数，配合 $β$ 的 $L_{1}$ 或 $L_{\infty}$ 惩罚项。另外一个有趣的例子，支持向量机中使用了“合页（hinge）”损失函数，这时的损失函数为分段线性函数而惩罚项为二次函数，但却产生了对偶空间（dual space）上分段线性的解路径，在第 12.3.5 节会更详细讨论。

3.8.3 丹齐格选择器

Candes and Tao (2007) 提出了下面的最小化问题：

\begin{matrix} (3.78) & min_{β} ‖ β ‖_{1} subject to ‖ X^{T} (y - X β) ‖_{\infty} \leq s \end{matrix}

并将其解称为 丹齐格选择器（Dantzig selector, DS）。可等价地写为：

\begin{matrix} (3.79) & min_{β} ‖ X^{T} (y - X β) ‖_{\infty} subject to ‖ β ‖_{1} \leq t \end{matrix}

其中的 $| \cdot |_{\infty}$ 为 $L_{\infty}$ 范式，即这个向量中的元素的最大绝对值。与套索回归的定义相比，这里用梯度的最大绝对值，代替了误差平方和本身。注意随着 t 的增大，若 $p < N$ ，两者均会得到最小二乘解；若 $p \geq N$ ，则两者都都会得到最小 $L_{1}$ 范式的最小二乘解。然而当 t 的取值比较小时，丹齐格选择器的解与套索回归不再一致。

Candes and Tao (2007) 证明了丹齐格选择器的解是一个线性规划（linear programming）问题，这也是其命名的来源，线性规划的单纯形法（simplex method）的发明者乔治丹齐格（George Dantzig）。论文中证明了这个方法很多特别的数学性质，使其可以计算出稀疏的隐含系数向量。Bickel et al. (2008) 证明了套索回归也有同样的性质。

但丹齐格选择器在实际操作中的性质却不尽如人意。其在思路上似乎与套索回归相似，尤其是在对照于套索回归的稳态条件 3.58。套索回归，与最小角回归类似，每一步中所有选入模型的变量与当前的残差有相同的内积（或相关程度），再将它们的系数向最小化残差平方和的方向移动。在整个移动过程中，其相同的内积（相关程度）单调地下降（练习 3.23），且这个值都大于未被选入模型的变量与残差的内积（相关程度）。丹齐格选择器与之不同，试图对当前残差与所有²自变量的内积的最大值进行最小化。因此其得到的最小化的最大值，要小于套索回归中对应的值，却也产生了和独特的解的性质：假设被选入模型的变量个数为 m，则在所有变量中会有 m 个变量与残差的相关程度等于解中的最大值，但是这 m 个变量却不一定都被选入模型。或者可以说，可能存在某个被选入模型的变量，其与当前残差的相关程度小于某个没有被选入模型的变量（Efron et al., 2007）。这样的解不太合常理，有时可能会使影响最终模型预测的准确度。Efron et al. (2007)也证明了丹齐格选择器可能产生极不规则的随正则化参数 $s$ 变化的系数路径。

3.8.4 分组套索回归

在某些场景中，部分自变量可能自然地有分组属性。例如，属于同一个生物途径（biological pathway）中的基因，或一个多类型的分类变量产生的一系列哑变量。这时可能希望对这一组中所有变量的系数进行同步的收缩。分组套索回归（grouped lasso） 可实现这一点。假设 $p$ 个自变量分属于 $L$ 个组，某个组 $ℓ$ 中包含的变量个数为 $p_{ℓ}$ 。为了书写简便，第 $ℓ$ 组的输入变量组成的矩阵写为 $X_{ℓ}$ ，对应的系数向量为 $β_{ℓ}$ 。则分组套索回归的最小化（凸函数）准则为：

\begin{matrix} (3.80) & min_{β \in R^{p}} (‖ y - β_{0} 1 - \sum_{ℓ = 1}^{L} X_{ℓ} β_{ℓ} ‖_{2}^{2} + λ \sum_{ℓ = 1}^{L} \sqrt{p_{ℓ}} ‖ β ‖_{2}) \end{matrix}

其中 $p_{ℓ}$ 为不同组中的变量个数， $‖ \cdot ‖_{2}$ 为（不加平方的）欧氏范式。由于当且仅当向量 $β_{ℓ}$ 中每个元素均为 0 时，其欧氏范式为 0，这个算法会倾向于拉大系数解在变量组间的距离，即对于某个 $λ$ 值，可能会使某个组中所有的变量都被排除出模型。Bakin (1999)、Lin and Zhang (2006) 提出了这个方法，Yuan 和 Lin（2007）对其继续进行研究和推广。其推广包括使用更广义的 $L_{2}$ 范式 $| η |_{K} = (η^{T} K η)^{1 / 2}$ ，或允许有交叉的变量分组（Zhao et al., 2008）。此方法可以与稀疏加性模型的拟合方法相关联（Lin and Zhang, 2006; Ravikumar et al., 2008）。

3.8.5 套索回归的其他性质

有很多对 $N$ 和 $p$ 很大时套索回归和其他相关算法对真实模型的拟合程度的研究。例如：Knight and Fu (2000)，Greenshtein and Ritov (2004)，Tropp (2004)，Donoho (2006b)，Meinshausen (2007)，Meinshausen and Bühlmann (2006)，Tropp (2006)，Zhao and Yu (2006)，Wainwright (2006)，Bunea et al. (2007)。Donoho (2006b) 研究了当 $p > N$ 时随着约束条件 $t$ 的增大套索回归解的变化。趋于极限的情况下的解是所有使得模型在训练集上拟合误差为 0 的解中 $L_{1}$ 范式最小的那个。他证明了在模型输入变量矩阵 $X$ 的特定假设下，如果真实模型的比较稀疏（sparse）³，则其解很可能会选择正确的输入变量进入模型。

这个领域中很多研究都对模型的输入变量矩阵有下面这种类型的假设：

\begin{matrix} (3.81) & max_{j \in S^{c}} ‖ x_{j}^{T} X_{S} (X_{S}^{T} X_{S})^{- 1} ‖_{1} \leq (1 - ϵ), for some ϵ \in (0, 1] \end{matrix}

其中 $S$ 为潜在的真实模型中的系数非零的输入特征变量集合， $X_{S}$ 为这些真实变量组成的输入变量矩阵； $S^{c}$ 为潜在的真实模型中的系数为零的输入特征变量集合， $X_{S^{c}}$ 则为相应的输入变量矩阵。此假设的含义为 $X_{S^{c}}$ 对 $X_{S}$ 的最小二乘回归系数足够小，即“有效”的输入变量 $S$ 与无效的变量 $S^{c}$ 相关程度足够小。

对于选入模型的变量，套索回归的收缩性会使得这些变量的系数估计值向 0 的方向偏差，而且通常一致⁴的估计。一个减少这个偏差的方法是先用套索回归识别出含有非零系数的变量，再用选中的这些特征变量进行无约束条件的线性回归。如果选中的变量过多，这个方法有时会不可行。另一种想法是嵌套地使用套索回归：用套索回归选择非零系数的自变量后，再用被选中的这些变量进行第二次套索回归。这即为（Meinshausen, 2007）提出的 松弛套索（relaxed lasso） 回归。两步套索回归的惩罚项参数均由交叉验证来决定。第一步的套索回归排除掉一部分无效的“干扰”变量，所以第二步的交叉验证常常会选择一个比较小的参数 $λ$ ，也因此对系数的收缩程度要低于第一步。

或者，可以对套索回归中的惩罚项进行修改，降低其对大系数的收缩程度。例如 Fan and Li (2005) 提出的 平滑裁剪绝对偏差（smoothly clipped absolute deviation，SCAD） 用 $J_{a} (β, λ)$ 代替 $λ | β |$ ，使得：

\begin{matrix} (3.82) & \frac{d J_{a} (β, λ)}{d β} = λ \cdot sign (β) [I (| β | \leq λ) + \frac{(a λ - | β |)_{+}}{(a - 1) λ} I (| β | > λ)] \end{matrix}

其中 $a \geq 2$ 。方括号中的第二项起到了减少套索回归对大系数值的收缩程度的作用，当 $a \to \infty$ 时，完全没有收缩。图 3.20 对比了不同的惩罚项函数：SCAD、套索回归、和 $| β |^{1 - v}$ 。然而 SCAD 的一个缺点是其惩罚项是非凸函数，导致其更难求解。（Zou, 2006）提出的 自适应套索（adaptive lasso） 回归使用一个加权的惩罚项 $\sum_{j = 1}^{p} w_{j} | β_{j} |$ ，其权重为 $w_{j} = 1 / | {\hat{β}}_{j} |^{v}$ ， ${\hat{β}}_{j}$ 为普通最小二乘的估计值，参数 $v > 0$ 。这实际上是对第 3.4.3 节中的 $| β |^{q}$ 惩罚项的近似（ $q = 1 - v$ ）。自适应套索回归的结果是一致估计量，而且其最优化准则函数与套索回归一样是凸函数。

**图 3.20**：套索回归惩罚项函数，和其他两种降低对大系数收缩程度的方法的非凸惩罚项函数。SCAD 的惩罚项函数系数为 $\lambda=1$ 和 $a=4$，最右图中 $v=\frac{1}{2}$。 — **图 3.20**：套索回归惩罚项函数，和其他两种降低对大系数收缩程度的方法的非凸惩罚项函数。SCAD 的惩罚项函数系数为 $λ = 1$ 和 $a = 4$ ，最右图中 $v = \frac{1}{2}$ 。

3.8.6 逐路径坐标优化

除了最小角回归（LAR）算法⁵外，另一种套索回归求解的方法为简单的 坐标下降（coordinate descent）。Fu (1998) 和 Daubechies et al. (2004) 提出了这个想法，Friedman et al. (2007) 和 Wu and Lange (2008) 以及其他作者进行了后续的研究和推广。其做法是在拉格朗日定义式中，给定惩罚项参数 $λ$ ，依次对单个的参数求解最优（而将其他参数值固定在当前步骤的取值上）。

假设自变量均已标准化为均值 0 和范式（方差）1。惩罚项系数为 $λ$ ，当前步骤中的 $β_{k}$ 估计值为 ${\tilde{β}}_{k} (λ)$ ，则可以将 $β_{j}$ 在定义 3.52 的准则函数中独立出来：

\begin{matrix} (3.83) & \begin{aligned} R (\tilde{β} (λ), β_{j}) = & \frac{1}{2} \sum_{i = 1}^{N} {(y_{i} - \sum_{k \neq j} x_{i k} {\tilde{β}}_{k} (λ) - x_{i j} β_{j})}^{2} + \\ λ \sum_{k \neq j} | {\tilde{β}}_{k} (λ) | + λ | β_{j} | \end{aligned} \end{matrix}

此处省略了截距项，其中的因子 $\frac{1}{2}$ 是为了后面推导的方便。上式可视为单变量的套索回归问题，其输出变量为偏残差（partial residual） $y_{i} - {\tilde{y}}_{i}^{(j)} = y_{i} - \sum_{k \neq j} x_{i k} {\tilde{β}}_{k} (λ)$ 。这个问题可以找到显式解，于是得到了这样的系数更新：

\begin{matrix} (3.84) & {\tilde{β}}_{j} (λ) \leftarrow S (\sum_{i = 1}^{N} x_{i j} (y_{i} - {\tilde{y}}_{i}^{(j)}), λ) \end{matrix}

其中的 $S (t, λ) = sign (t) (| t | - λ)_{+}$ 即为在表 3.4 中提到的软阈值算子，其中的第一个输入参数为偏残差对标准化的变量 $x_{i j}$ 的单变量最小二乘系数估计。对所有变量依次进行这个更新，并循环这个迭代过程，直到满足某个收敛条件，最终可得到套索回归的估计 $\hat{β} (λ)$ 。

这个简单的算法可以在 $λ$ 的取值网格上快速地求解套索回归。先从最大⁶的参数值 $λ_{max}$ 开始，此时的系数估计为 $\hat{β} (λ_{max}) = 0$ 。然后逐渐降低 $λ$ 的值，对每个 $λ$ ，循环迭代系数的更新过程，直到满足收敛条件。在上述的每一步中，将之前一步得到的解作为“热启动”的初始值，即从上一步 $β$ 的解开始更新过程。这种算法，尤其在高维度问题中，可能要快于最小角回归算法。其计算效率归功于可以快速计算不同 j 的等式 3.84 更新，以及 ${\tilde{β}}_{j}$ 会在很多步骤中保持为 0。但另一方面，这个算法得到的是 $λ$ 取值在某个离散的网格上的解，而不是完整的解的路径（ $β$ 对 $λ$ 的函数）。这种算法也可以应用在弹性网络（elastic net）、分组套索回归、和其他惩罚项可分解为各个系数的惩罚项加总的模型上（Friedman et al., 2010）。通过大量修正后，这个算法也可应用在融合套索（fused lasso）回归上（第 18.4.2 节）；更多细节见 Friedman et al. (2007)。

原文脚注 4：在第一版中，该算法为第十章的算法 10.4。 ↩︎
即包括当前选入模型的以及没有选入模型的所有输入变量。 ↩︎
译者不确定如何理解“if the true model is sparse”，可能指的是真实模型的不同变量的系数值差异比较大。 ↩︎
原文脚注 5：统计学的“一致性”是指当随着样本量增大，估计值收敛到真实值。 ↩︎
原文为“LARS”，这是第一次出现这个缩写，在后面章节中也有出现。译者认为可能不同的作者采取了不同的缩写方式，这里指的是最小角回归（LAR）。 ↩︎
原文为“最小”，这与后面脚标中的 max 矛盾，译者理解此处应该为“最大”，对应着最强的约束条件。 ↩︎