13 稀疏线性模型¶

13.1 引言¶

我们在3.5.4节中介绍了关于特征选择的主题，讨论了寻找与输出具有高互信息的输入变量的方法。这种方法的问题在于，它基于一种短视的策略，即一次只看一个变量。如果输入变量之间存在交互作用，则此操作可能会失败。例如，如果$y = {\rm{xor}}(x_1, x_2)$,那么$x_1$和$x_2$都不能单独预测响应，但是它们一起工作就可以进行完美地预测。举一个真实的例子，考虑一下基因关联研究:有时两个基因本身可能是无害的，但当它们一起出现时，就会导致隐性疾病。

在本章中，我们将使用基于模型的方法一次选择一组变量。如果模型是广义线性模型，对于某些形式为$p(y|\mathbf{x})=p(y|f(\mathbf{w^Tx}))$的连接函数$f$，则可以通过使权重向量$\mathbf{w}$稀疏来进行特征选择，即，向量中很多元素为0。这种方法提供了重要的计算优势，这一点我们将在下文展示。

以下是一些特征选择/稀疏化有用的应用：

在许多问题中，我们有比训练样本数量$N$更多的特征维度$D$。在这种情况下，设计矩阵会显得又矮又胖，而不是又高又苗条。这被称为小$N$大$D$问题。当我们开发更多的高吞吐量测量设备时，这种情况变得越来越普遍，例如，使用基因微阵列，通常测量$D∼10,000$个基因的表达水平，但是只得到$N∼100$个这样的样例。(这也许是一个时代的迹象，甚至我们的数据似乎也在变胖…)。我们可能想要找到能够准确预测响应的最小特征集(例如，细胞的生长速度)，以防止过度拟合，降低制造诊断设备的成本，或帮助科学地洞察问题。
在第14章，我们将使用以训练样本为中心的基函数，即$\phi(\mathbf{x})=[\mathcal{k}(\mathbf{x},\mathbf{x}_1),...,\mathcal{k}(\mathbf{x},\mathbf{x}_N)]$，其中$\mathcal{k}$为一个核函数。通过这种方式得到的设计矩阵大小为$N\times N$。在这种情况下的特征选择等价于选择训练样本的一个子集，从而缓解过拟合和减少计算成本。这被称为稀疏核机器(sparse kernel machine)。
在信号处理中，通常会使用小波基函数来表示信号(图像，声音等等)。为了节约时间和空间，寻找关于信号的稀疏表示是十分有用的，也就是使用少量的基函数。这将允许我们使用少量的测量装置来实现对信号的估计，也就是对信号进行了压缩。13.8.3将介绍更多的信息。

需要注意的是，特征选择和稀疏化是目前机器学习/统计领域最活跃的领域之一。本章，我们将只给出一些主要结论的概述。

13.2 贝叶斯观点下的变量选择¶

接下来，我们将介绍解决变量选择问题的一个很自然的方法。如果特征$j$被认为是“相关的”，那么令$\gamma_j=1$，否则令$\gamma_j=0$。我们的目标是计算模型的后验分布 $$ p(\mathbf{\gamma}|\mathcal{D})=\frac{e^{-f(\gamma)}}{\sum_{\gamma^\prime}e^{-f(\gamma^\prime)}} \tag{13.1} $$ 其中$f(\gamma)$是成本函数： $$ f(\gamma)\triangleq -[\log p(\mathcal{D|\gamma})+\log p(\gamma)] \tag{13.2} $$ 举例来说，假设我们从$D=10$的线性回归模型$y_i \sim \mathcal{N}(\mathbf{w}^T\mathbf{x}_i,\sigma^2)$中产生了$N=20$个样本，其中权重向量$\mathbf{w}$中的$K=5$个元素为非零元素。特别地，我们使用$\mathbf{w}=(0.00,-1.67,0.13,0.00,0.00,1.19,0.00,-0.04,0.33,0.00)$以及$\sigma^2=1$。我们可以枚举出所有的$2^{10}=1024$个模型，并且对每个模型计算$p(\gamma|\mathcal{D})$。（我们将在后文介绍如何计算）。我们使用**格雷码(Gray code)**的顺序对所有模型进行排序，从而保证连续的模型之间只相差一个比特位（这么做的理由是考虑计算上的优势，我们将在13.2.3节讨论这一点）。

图13.1(a)展示了最终的比特模式图。每个模型的成本$f(\gamma)$如图13.1(b)所示。我们发现目标函数极度地“不平衡”。如果我们计算出关于模型的后验分布$p(\gamma|\mathcal{D})$，那么这种情况就可以很容易地解释了。后验分布如图13.1(c)所示。排名前8的模型如下表所示：

model	prob	members
4	0.447	2,
61	0.241	2,6,
452	0.103	2,6,9,
60	0.091	2,3,6,
29	0.041	2,5,
68	0.021	2,6,7,
36	0.015	2,5,6,
5	0.010	2,3,

“真正”的模型应该是$\{2,3,6,8,9\}$。然而，与特征$3$和特征$8$关联的系数却十分的小（相对于$\sigma ^ 2$），所以这些变量很难被检测到。给定足够的数据，这种方法将收敛于真正的模型（假设数据产生于同一个线性模型），但对于有限的数据集而言，后验分布通常存在非常大的不确定性。

当后验分布所基于的模型数量十分大时，解释起来会十分困难，所以我们寻找各种具有总结意义的统计量。一个自然的选择是后验分布的峰值，即MAP估计： $$ \hat{\mathbf{\gamma}}=\arg \max p(\gamma|\mathcal{D})=\arg \min f(\gamma) \tag{13.3} $$ 然而，对于整个后验分布的概率质量而言，峰值通常不具备代表性（见5.2.1.3节）。一种更好的总结方法是**中值模型(median model)**，即： $$ \hat{\mathbf{\gamma}}=\{j:p(\gamma_j=1|\mathcal{D})>0.5\} \tag{13.4} $$ 这要求我们计算后验边缘**包含概率(inclusion probabilities)** $p(\gamma_j=1|\mathcal{D})$。这些结果在图13.1(d)中进行展示。我们发现如果变量2和6被包含，那么模型会很有信心；如果我们将决策阈值降低到0.1，模型将会同时考虑变量3和9.然而，如果我们还想“获取”变量8，那么于此同时会将2个错误的正例(5和7)囊括进来。这种在假正例和假负例之间进行权衡的细节在5.7.2.1节中已有讨论。

上述例子说明了一个变量选择的“金本体制(gold standard)”：当所涉及的问题足够小（只有10个变量）时，我们可以精确地计算出完整的后验分布。当然，变量选择对于维度数量特别大的情况是最有用的。因为存在$2^D$个可能的模型（位向量），通常情况下不可能计算出完整的后验分布，以及一些统计量，比如MAP估计或者边缘包含概率。因此，本章我们将集中介绍如何进行算法的加速。但在此之前，我们将介绍在上面的例子中，是如何计算$p(\gamma|\mathcal{D})$。

13.2.1 spike and slab model¶

后验分布定义为： $$ p(\mathbf{\gamma}|\mathcal{D}) \propto p(\mathbf{\gamma})p(\mathcal{D}|\mathbf{\gamma}) \tag{13.5} $$ 我们首先考虑先验分布，其次是似然函数。

通常情况下，我们使用如下的关于位向量的先验分布： $$ p(\mathbf{\gamma})=\prod_{j=1}^D{\rm{Ber}}(\gamma_j|\pi_0)=\pi_0^{\parallel\mathbf{\gamma}\parallel_0}(1-\pi_0)^{D-\parallel\mathbf{\gamma}\parallel_0} \tag{13.6} $$ 其中$\pi_0$表示一个特征是相关联的概率，$\parallel\mathbf{\gamma}\parallel_0=\sum_{j=1}^D\gamma_j$为$\mathcal{l}_0$**伪范数(pseudo-norm)**，即向量中非零元素的数量。为了方便与后面模型进行比较，对上式取对数形式是十分有用的： $$

(14)¶\[\begin{align} \log p(\mathbf{\gamma}|\pi_0) & = \parallel\mathbf{\gamma}\parallel_0 \log \pi_0+(D-\parallel\mathbf{\gamma}\parallel_0)\log(1-\pi_0) \tag{13.7} \\ & =\parallel\mathbf{\gamma}\parallel_0 (\log \pi_0-\log(1-\pi_0))+{\rm{const}} \tag{13.8}\\ & = - \lambda \parallel\mathbf{\gamma}\parallel_0 + {\rm{const}} \tag{13.9} \end{align}\]

\[ \begin{align}\begin{aligned} 其中$\lambda \triangleq \log \frac{1-\pi_0}{\pi_0} $控制着模型的稀疏强度。\\我们可以按照下面的方式定义似然函数： \end{aligned}\end{align} \]

p(\mathcal{D}|\mathbf{\gamma})=p(\mathbf{y}|\mathbf{X,\gamma})=\iint p(\mathbf{y}|\mathbf{X,w,\gamma})p(\mathbf{w}|\mathbf{\gamma},\sigma^2)p(\sigma^2)d\mathbf{w}d\sigma^2 \tag{13.10} $$ 出于符号上的方便，我们假设响应值已经被中心化处理，（即$\bar{y}=0$）,所以我们可以忽略所有的偏置项$\mu$。

我们现在讨论先验分布$p(\mathbf{w}|\mathbf{\gamma},\sigma^2)$。如果$\gamma_j=0$，即特征$j$是无关的，此时我们期望$w_j=0$。如果$\gamma_j=1$,我们期望$w_j$是非零的值。如果我们对输入进行标准化处理，一个合理的先验分布应该是$\mathcal{N}(\mathbf{0},\sigma^2\sigma_w^2)$，其中$\sigma_w^2$控制着我们期望系数与相关性变量（其变化范围由全局噪音水平$\sigma^2$表征）的相关程度。我们可以将这个先验分布总结为： $$ p(w_j|\sigma^2,\gamma_j)= \begin{cases} \delta_0(w_j), & \mbox{if } \gamma_j=0 \\ \mathcal{N}(w_j|0,\sigma^2\sigma_w^2) & \mbox{if }\gamma_j=1 \tag{13.11} \end{cases} $$ 第一项为在原点处的一个"spike(针)"。当$\sigma_w^2 \rightarrow \infty$，分布$p(w_j|\gamma_j=1)$将趋向于一个均匀分布，这被认为是一个具有恒定高度的“slab（平板）”。所以这被称为spike and slab模型。

我们可以将那些$w_j=0$的系数从模型中“清理”处理，因为在先验分布中，它们的取值被限制为0。所以式13.10变成如下的形式（假设似然函数服从高斯分布）： $$ p(\mathcal{D}|\mathbf{\gamma})=\iint \mathcal{N}(\mathbf{y}|\mathbf{X}_\gamma\mathbf{w}_\gamma,\sigma^2\mathbf{I}_N)\mathcal{N}(\mathbf{w}_\gamma|\mathbf{0}_\gamma,\sigma^2\sigma_w^2\mathbf{I}_\gamma)p(\sigma^2)d\mathbf{w}_\gamma d\sigma^2 \tag{13.12} $$ 其中$\mathbf{X}\gamma$为新的设计矩阵，其中我们只选择了那些$\gamma_j=1$的列的值，$\mathbf{0}\gamma=\mathbf{0}{\parallel\gamma\parallel_0},\mathbf{I}\gamma=\mathbf{I}{\parallel\gamma\parallel_0}$。类似地，我们通过符号$\mathbf{w}\gamma$来提醒读者该系数只表示那些$\gamma_j=1$的情况。在后面的内容中，我们定义一种更具一般性的先验分布$p(\mathbf{w}|\gamma,\sigma^2)=\mathcal{N}(\mathbf{w}|\mathbf{0}\gamma,\sigma^2\mathbf{\Sigma}\gamma)$，其中$\mathbf{\Sigma}\gamma$表示任意的正定矩阵，而无需要求它满足$\mathbf{\Sigma}\gamma=\sigma_w^2\mathbf{I}_\mathbf{\gamma}$.

基于上述的先验分布，我们现在可以计算边缘似然。如果噪声的方差已知，我们可以计算出（使用公式13.151）如下的边缘似然： $$

(15)¶\[\begin{align} p(\mathcal{D}|\gamma,\sigma^2) & = \int\mathcal{N}(\mathbf{y}|\mathbf{X}_\gamma\mathbf{w}_\gamma,\sigma^2\mathbf{I})\mathcal{N}(\mathbf{0},\sigma^2\mathbf{\Sigma}_\gamma)d\mathbf{w}_\gamma=\mathcal{N}(\mathbf{y}|\mathbf{0},\mathbf{C}_\gamma) \tag{13.13}\\ \mathbf{C}_\gamma & \triangleq \sigma^2\mathbf{X}_\gamma\mathbf{\Sigma}_\gamma\mathbf{X}_\gamma^T+\sigma^2\mathbf{I}_N \tag{13.14} \end{align}\]

\[ 如果噪声是未知的，我们可以为其添加一个先验分布，并且对它进行积分。通常我们使用先验分布$p(\sigma^2)={\rm{IG}}(\sigma^2|a_\sigma,b_\sigma)$。关于超参$a,b$的设置在一些文献中可以找到。如果我们使用$a=b=0$，我们将使用Jeffrey先验分布，即$p(\sigma^2)\propto \sigma^{-2}$。当我们对噪音项进行积分，我们将得到如下的更加复杂的关于边缘似然的表达式： \]

\[ \begin{align}\begin{aligned}## 13.3 $\mathcal{l}_1$ 正则：基本原理\\当我们的变量数量特别多时，计算后验分布$p(\gamma|\mathcal{D})$的峰值会十分困难。尽管贪心算法通常会奏效，但同时会陷入局部最优解。\\部分原因是因为变量$\gamma_j$是离散的，即$\gamma_j \in \{0,1\}$。在最优化领域，通常会放松这种hard级别的约束，使用连续变量代替离散变量。在前文介绍的spike-and-slab形式的先验分布中，事件$w_j=0$被赋予的概率质量是有限的，我们可以使用一种连续的先验分布替代它，在这种先验分布中，我们通过将更多的概率密度赋予0周边的区域，从而“鼓励”$w_j=0$。这种想法在7.4节的鲁棒性线性回归中首次被介绍。在那里我们介绍了Laplace分布具有肥尾的事实。此处，我们将探索该分布在$\mu=0$的附近呈现针状的事实。更加精确的表达形式为： \end{aligned}\end{align} \]

p(\mathbf{w}|\lambda)=\prod_{j=1}^D{\rm{Lap}}(w_j|0,1/\lambda) \propto \prod_{j=1}^De^{-\lambda|w_j|} \tag{13.32} $$ 我们将赋予偏置量$w_0$一个均匀分布，即$p(w_0) \propto 1$。我们基于上述先验分布求解MAP估计。含惩罚项的负对数似然形式为： $$ f(\mathbf{w})=-\log p(\mathcal{D}|\mathbf{w})-\log p(\mathbf{w}|\lambda)={\rm{NLL}}(\mathbf{w})+\lambda||\mathbf{w}||_1 \tag{13.33} $$ 其中$||\mathbf{w}||1=\sum{j=1}^D|w_j|$表示$\mathbf{w}$的$\mathcal{l}_1$正则。当$\lambda$的值比较合适时，估计量$\hat{\mathbf{w}}$将会变得稀疏，其原因我们将会在后文介绍。的确这可以看作是非凸$\mathcal{l}_0$目标的凸近似 $$ \mathop{\arg\min}_{\mathbf{w}} {\rm{NLL}}(\mathbf{w})+\lambda||\mathbf{w}||_0 \tag{13.34} $$ 在线性回归中，$\mathcal{l_1}$目标函数为： $$

(16)¶\[\begin{align} f(\mathbf{w}) & = \sum_{i=1}^N -\frac{1}{2\sigma^2}(y_i-(w_0+\mathbf{w}^T\mathbf{x}_i))^2+\lambda||\mathbf{w}||_1 \tag{13.35} \\ & = {\rm{RSS}}(\mathbf{w})+\lambda^\prime||\mathbf{w}||_1 \tag{13.36} \end{align}\]

\[ \begin{align}\begin{aligned} 其中$\lambda^\prime=2\lambda\sigma^2$。这种方法被称为**基追踪去噪(basis pursuit denoising, BPDN)**。之所以叫这个名字我们将在后面介绍。一般情况下，为参数赋予一个均值为0的Laplace先验分布并计算MAP估计的技术被称为$\mathcal{l}_1$**正则(regularization)**。它可以与任意的凸或者非凸NLL项结合。许多针对该问题的算法被设计出来，我们将在13.4节介绍其中的一些内容。\\### 13.3.1 为什么$\mathcal{l}_1$正则可以得到稀疏解\\现在我们来解释为什么$\mathcal{l}_1$正则可以导致一个稀疏的解，而$\mathcal{l}_2$不行。我们以线性回归为例，尽管相似的结论同样适用于逻辑回归和其他广义线性模型。\\BPDN目标函数是具备如下形式的非平滑目标函数： \end{aligned}\end{align} \]

\mathop{\min}_{\mathbf{w}} {\rm{RSS}}(\mathbf{w})+\lambda||\mathbf{w}||_1 \tag{13.37} $$ 我们可以将上式改写成含约束的光滑的目标函数形式： $$ \mathop{\min}_{\mathbf{w}} {\rm{RSS}}\ s.t. ||\mathbf{w}||_1 \le B \tag{13.38} $$ 其中$B$为权重的$\mathcal{l}_1$范数的上确界：一个小的（紧的）边界值$B$对应于一个大的惩罚系数$\lambda$，反之亦然。式13.38被称为lasso，即”最小绝对收缩和选择操作(least absolute shrinkage and selection operator)”。我们将在后文介绍这个名字的由来。

类似地，我们可以将岭回归的形式写成： $$ \mathop{\min}_{\mathbf{w}} {\rm{RSS}}(\mathbf{w})+\lambda||\mathbf{w}||_2^2 \tag{13.39} $$ 或者一种含边界约束的版本： $$ \mathop{\min}_{\mathbf{w}} {\rm{RSS}}(\mathbf{w}) s.t. ||\mathbf{w}||_2^2 \le B \tag{13.40} $$ 在图13.3中，我们绘制了$\rm{RSS}$目标函数以及$\mathcal{l}_2$和$\mathcal{l}_1$约束表面的轮廓线。根据含约束优化的理论，我们知道当目标函数的最低集合与约束表面（假设约束是激活的）相交时，才会出现最优解。如果我们放松约束$B$，那么这一点在几何上会显得更加清晰，我们“放大”$\mathcal{l}_1$“球”直到它与目标接触；显然相较于“边”，“球”的角要更容易与椭圆接触，尤其是高维的情况，因为角更加突出。角对应着稀疏解，它坐落在坐标轴上。相反，当我们放大$\mathcal{l}_2$球的时候，它可以在任何点与目标接触；因为没有角的存在，所以也没有稀疏化的倾向。

为了从另一个角度说明这一点，当我们使用岭回归时，稀疏解，比如$\mathbf{w}=(1,0)$的先验成本，与稠密解，比如$\mathbf{w}=(1/\sqrt{2},1/\sqrt{2})$，只要它们有同样的$\mathcal{l}_2范数$： $$ ||(1,0)||_2=||(1/\sqrt{2},1/\sqrt{2})||_2=1 \tag{13.41} $$ 然而，对于lasso,$\mathbf{w}=(1,0)$要比$\mathbf{w}=(1/\sqrt{2},1/\sqrt{2})$便宜，因为： $$ ||(1,0)||_1=1\lt||(1/\sqrt{2},1/\sqrt{2})||_1=\sqrt{2} \tag{13.42} $$ 说明$\mathcal{l}_1$正则能导致稀疏解的检验最优解成立时的条件。我们将在13.3.2节介绍。 $$ Figure 13.4 $$

13.3.2 lasso最优解条件¶

lasso的目标函数形式为： $$ f(\mathbf{\theta})={\rm{RSS}}(\mathbf{\theta})+\lambda||\mathbf{w}||_1 \tag{13.43} $$ 不幸地是，当$w_j=0$时，$||\mathbf{w}||_1$项是不可微的。这是一个**非平滑(non-smooth)**优化问题的例子。

为了解决非平滑函数，我们需要拓展导数的符号。我们定义一个（凸）函数$f:\mathcal{I}\rightarrow \mathbb{R}$在点$\theta_0$的**次导数(subderivative)或者次梯度(subgradient)**为一个标量$g$，使其满足: $$ f(\theta)-f(\theta_0)\ge g(\theta-\theta_0) \ \forall\theta\in\mathcal{I} \tag{13.44} $$

其中$\mathcal{I}$为包含$\theta_0$的区间。图13.4给出了示意图。我们定义次导数的集合为区间$[a,b]$，其中$a,b$分别为单侧极限 $$ a=\mathop{\lim}_{\theta\rightarrow\theta_0^-}\frac{f(\theta)-f(\theta_0)}{\theta-\theta_0},b=\mathop{\lim}_{\theta\rightarrow\theta_0^+}\frac{f(\theta)-f(\theta_0)}{\theta-\theta_0} \tag{13.46} $$ 所有次导数的集合$[a,b]$被称为函数$f$在点$\theta_0$处的**次微分(subdifferential)**，符号上表示为$\partial f(\theta)|{\theta_0}$。举例来说，绝对值函数$f(\theta)=|\theta|$，其次导数定义为： $$ \partial f(\theta)= \begin{cases} \{-1\}, & {\rm{if}}\ \theta\lt0 \\ [-1,1], & {\rm{if}}\ \theta=0 \\ \{+1\}, & {\rm{if}}\ \theta\gt0 \tag{13.47} \end{cases} $$ 如果函数在每个点都可微，则$\partial f(\theta)={\frac{df(\theta)}{d\theta}}$。类别于标准的计算结果，函数$f$局部最小值所在的点$\hat{\theta}$满足充要条件$0\in\partial f(\theta)|\hat{\theta}$。

接下来我们将上述概念应用在lasso问题中。首先忽略非平滑的惩罚项。结果表明： $$

(17)¶\[\begin{align} \frac{\partial}{\partial w_j} &=a_jw_j-c_j \tag{13.48} \\ a_j & = 2\sum_{i=1}^{n}x_{ij}^2 \tag{13.49} \\ c_j & = 2\sum_{i=1}^{n}x_{ij}(y_i-\mathbf{w}_{-j}^T\mathbf{x}_{i,-j}) \tag{13.50} \end{align}\]

\[ \begin{align}\begin{aligned} 其中$\mathbf{w}_{-j}$表示除了分量$j$以外的元素，$\mathbf{x}_{i,-j}$表示相同的意思。我们发现$c_j$是（正比于）第$j$个特征$\mathbf{x}_{:,j}$与因为其他特征所导致的残差项$\mathbf{r}_{-j}=\mathbf{y}-\mathbf{X}_{:,-j}\mathbf{w}_{-j}$之间的相关系数。所以$c_j$的幅值表征了特征$j$对于预测$\mathbf{y}$的重要性（相对于其他特征和当前参数）。\\在增加了惩罚项之后，我们发现次导数由下式给定： \end{aligned}\end{align} \]

(18)¶\[\begin{align} \partial_{w_j}f(\mathbf{w}) & = (a_jw_j-c_j)+\lambda\partial_{w_j}||\mathbf{w}||_1 \tag{13.51} \\ & = \begin{cases} \{a_jw_j-c_j-\lambda\} & {\rm{if}}\ w_j \lt 0 \\ [-c_j-\lambda,-c_j+\lambda] & {\rm{if}} \ w_j=0 \\ \{a_jw_j-c_j+\lambda\} & {\rm{if}} \ w_j \gt 0 \tag{13.52} \end{cases} \end{align}\]

\[ 因此$\mathbf{w}$是局部最优解的充要条件是: \]

\mathbf{X}^T(\mathbf{X}\mathbf{w-\mathbf{y}})_j \in \begin{cases} {-\lambda} & {\rm{if}} \ w_j \lt 0 \ [-\lambda, \lambda] & {\rm{if}} \ w_j = 0 \ {\lambda} & {\rm{if}} \ w_j \gt 0 \tag{13.53} \end{cases} $$ 译者注:个人以为在式(13.53)中在两侧的导数符号反了。

根据$c_j$的取值，$\partial_{w_j}f(\mathbf{w})=0$的解存在如下的三种情况：

如果$c_j \lt -\lambda$，也就是说特征与残差呈强烈的负相关关系，在这种情况下次导数为0的条件是$\hat{w}_j=\frac{c_j+\lambda}{a_j}\lt0$.
如果$c_j \in [-\lambda,\lambda]$，也就是说特征与残差之间呈现弱相关的关系，此时次导数为0的条件是$\hat{w}_j=0$.
如果$c_j \gt \lambda$,也就是说特征与残差呈现强烈的正相关关系，此时次导数为0的条件是$\hat{w}_j=\frac{c_j-\lambda}{a_j}\gt 0$.

总结下来，我们有：

\[\begin{split} \hat{w}_j(c_j)= \begin{cases} (c_j+\lambda)/a_j & {\rm{if}} \ c_j \lt -\lambda \\ 0 & {\rm{if}} \ c_j \in [-\lambda, \lambda] \\ (c_j - \lambda)/a_j & {\rm{if}}\ c_j \gt \lambda \tag{13.54} \end{cases} \end{split}\]

我们可以将上式写成: $$ \hat{w}_j = {\rm{soft}}(\frac{c_j}{a_j};\frac{\lambda}{a_j}) \tag{13.55} $$ 其中 $$ {\rm{soft}}(a;\delta)\triangleq {\rm{sign}}(a)(|a|-\delta)_+ \tag{13.56} $$ 其中$x_+=\max(x,0)$。上式被称为**柔性阈值(soft thresholding)**。图13.5(a)给出了说明,其中我们绘制了$\hat{w}_j$与$c_j$的关系曲线。黑色虚线$w_j=c_j/a_j$对应于最小二乘解。红色实线，代表了正则估计$\hat{w}_j(c_j)$，该解将虚线向下(或向上)移动了$\lambda$,除了当$-\lambda \le c_j \le \lambda$，此时的$\hat{w}_j=0$。

相反，在图13.5(b)中，我们展示了生硬阈值(hard thresholding)。在这种情况下，当$-\lambda \le c_j \le \lambda$时，$w_j=0$，但在这个区间之外的$w_j$并没有被压缩。柔性阈值的曲线并没有与对角线相一致，这意味着哪怕是很大的系数都有向0收缩的趋势；因此lasso是一个有偏估计。这并不是我们想要的，因为如果似然值（通过$c_j$）表明系数$w_j$应该很大，我们并不希望对这些参数进行收缩。我们将在13.6.2节讨论关于这个话题的更多细节。

最后我们可以理解为什么Tibshirani发明了单词”lasso”:它的全称是”least absolute selection and shrinkage operator”,因为它在所有变量中选择了一个子集，并且通过惩罚绝对值的方式对所有系数进行了收缩。如果$\lambda=0$，我们将得到$OLS$解。如果$\lambda \ge \lambda_{max}$，我们将得到$\hat{\mathbf{w}}=0$,其中 $$ \lambda_{max}=||\mathbf{X}^T\mathbf{y}||_{\infty}=\max_j|\mathbf{y}^T\mathbf{x}_{:,j}| \tag{13.57} $$ 上式的计算依据为：如果对于所有的$j$，$(\mathbf{X^T\mathbf{y}})_j \in [-\lambda, \lambda]$都成立，那么$\mathbf{0}$将会是最优解。一般情况下，对于一个$\mathcal{l}_1$正则目标函数，其最大的惩罚为: $$ \lambda_{max}=\max_j|\nabla_jNLL(\mathbf{0})| \tag{13.58} $$

13.3.3 最小二乘，lasso,ridge和子集选择的比较¶

通过对比lasso与最小二乘,含$\mathcal{l}_2$正则和$\mathcal{l}_0$正则的最小二乘的区别，我们可以更加深入地理解$\mathcal{l}_1$正则。为了简单起见，假设$\mathbf{X}$所有的特征都已经被正交化，即$\mathbf{X}^T\mathbf{X}=\mathbf{I}$。在这种情况下，$RSS$由下式给定: $$

(19)¶\[\begin{align} {\rm{RSS}}(\mathbf{w}) & = ||\mathbf{y}-\mathbf{Xw}||^2=\mathbf{y}^T\mathbf{y}+\mathbf{w}^T\mathbf{X}^T\mathbf{X}\mathbf{w}-2\mathbf{w}^T\mathbf{X}^T\mathbf{y} \tag{13.59} \\ & = {\rm{const}} + \sum_kw_k^2 -2 \sum_k\sum_iw_kx_{ik}y_i \tag{13.60} \end{align}\]

\[ \begin{align}\begin{aligned} 所以我们发现，上式分解为了很多项的和，其中每个维度代表一项。所以，我们可以写出$MAP$和$ML$估计的解析解，具体如下：\\- $\rm{\mathbf{MLE}}$:$OLS$解为: \end{aligned}\end{align} \]

\hat{w}k^{OLS}=\mathbf{x}{:,k}^T\mathbf{y} \tag{13.61} $$

机器学习：概率视角

13 稀疏线性模型

Contents