ESL3.4 学习笔记（含奇异值分解与拉格朗日乘子法内容岭回归，Lasso回归，最小角回归及三者对比分析）

3.4 收缩的方法

这是一篇有关《统计学习基础》，原书名The Elements of Statistical Learning的学习笔记，该书学习难度较高，有很棒的学者将其翻译成中文并放在自己的个人网站上，翻译质量非常高，本博客中有关翻译的内容都是出自该学者的网页，个人解读部分才是自己经过查阅资料和其他学者的学习笔记，结合个人理解总结成的原创内容。

原文	The Elements of Statistical Learning
翻译	szcf-weiya
时间	2018-08-21
解读	Hytn Chen
更新	2020-02-17

文章目录3.4 收缩的方法翻译原文岭回归个人解读拉格朗日乘子法奇异值分解Lasso讨论：子集的选择，岭回归，Lasso最小角回归个人解读LAR 和 Lasso 自由度公式个人解读
超级长文预警~ 翻译原文

通过保留一部分预测变量而丢弃剩余的变量，子集选择 (subset selection) 可得到一个可解释的、预测误差可能比全模型低的模型．然而，因为这是一个离散的过程（变量不是保留就是丢弃），所以经常表现为高方差，因此不会降低全模型的预测误差．而收缩方法 (shrinkage methods) 更加连续，因此不会受 高易变性 (high variability) 太大的影响．

岭回归

岭回归 (Ridge regression) 根据回归系数的大小加上惩罚因子对它们进行收缩．岭回归的系数使得带惩罚的残差平方和最小

β^ridge=arg⁡min⁡β{∑i=1N(yi−β0−∑j=1pxijβj)2+λ∑j=1pβj2}(3.41) \hat{\beta}^{ridge}=\underset{\beta}{\arg\min}\Big\{\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2+\lambda\sum\limits_{j=1}^p\beta_j^2\Big\}\tag{3.41} β^​ridge=βargmin​{i=1∑N​(yi​−β0​−j=1∑p​xij​βj​)2+λj=1∑p​βj2​}(3.41)

这里λ≥0\lambda\ge 0λ≥0是控制收缩程度的参数：λ\lambdaλ值越大，收缩的程度越大．每个系数都向零收缩．通过参数的平方和来惩罚的想法也用在了神经网络，也被称作权重衰减 (weight decay)（第 11 章）．

岭回归问题可以等价地写成

β^ridge=arg⁡min⁡β∑i=1N(yi−β0−∑j=1pxijβj)2subject to ∑j=1pβj2≤t(3.42) \begin{aligned} \hat{\beta}^{ridge}&=\underset{\beta}{\arg\min}\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2\\ & \text{subject to }\sum\limits_{j=1}^p\beta_j^2 \le t \end{aligned} \tag{3.42} β^​ridge​=βargmin​i=1∑N​(yi​−β0​−j=1∑p​xij​βj​)2subject to j=1∑p​βj2​≤t​(3.42)

上式用参数显式表达了对回归参数大小的约束．

!!! note “weiya 注：”
式 (3.41)(3.41)(3.41) 其实是对式 (3.42)(3.42)(3.42) 应用 Lagrange 乘子法得到的．

(3.41)(3.41)(3.41) 中的 λ\lambdaλ 和 (3.42)(3.42)(3.42) 中的 ttt 存在一一对应．当在线性回归模型中有许多相关变量，它们的系数可能很难确定且有高方差．某个变量的较大的正系数可以与相关性强的变量的差不多大的负系数相互抵消．通过对系数加入大小限制，如 (3.42)(3.42)(3.42)，这个问题能得以减轻．

!!! note “weiya 注：”
这里说的是，在没有对参数大小进行限制前，会存在一对相关性强的变量，它们系数取值符号相反，但绝对值差不多大，会大大增加方差，这也就是高方差的体现，但其实它们的合作用效果近似为 000，所以考虑引进对参数大小的惩罚．

对输入按比例进行缩放时，岭回归的解不相等，因此求解 (3.41)(3.41)(3.41) 前我们需要对输入进行标准化．另外，注意到惩罚项不包含截距 β0\beta_0β0​．对截距的惩罚会使得过程依赖于 Y\mathbf{Y}Y 的初始选择；也就是，对每个 yiy_iyi​ 加上常数 ccc 不是简单地导致预测值会偏离同样的量 ccc．可以证明（练习 3.5）经过对输入进行中心化（每个 xijx_{ij}xij​ 替换为 xij−xˉjx_{ij}-\bar x_jxij​−xˉj​）后，(3.41)(3.41)(3.41) 的解可以分成两部分．我们用 yˉ=1N∑1Nyi\bar y=\frac{1}{N}\sum_1^Ny_iyˉ​=N1​∑1N​yi​ 来估计 β0\beta_0β0​．剩余的参数利用中心化的 xijx_{ij}xij​ 通过无截距的岭回归来估计．今后我们假设中心化已经完成，则输入矩阵 X\mathbf XX 有 ppp（不是 p+1p+1p+1）列．

!!! info “weiya 注：Ex. 3.5”
已解答，详细证明过程见 Issue 95: Ex. 3.5

将 (3.41)(3.41)(3.41) 的准则写成矩阵形式

RSS(λ)=(y−Xβ)T(y−Xβ)+λβTβ(3.43) \rm{RSS}(\lambda)=(\mathbf{y}-\mathbf{X}\beta)^T(\mathbf{y}-\mathbf{X}\beta)+\lambda\beta^T\beta \tag{3.43} RSS(λ)=(y−Xβ)T(y−Xβ)+λβTβ(3.43)

可以简单地看出岭回归的解为

β^ridge=(XTX+λI)−1XTy(3.44) \hat{\beta}^{ridge}=(\mathbf{X^TX}+\lambda \mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.44} β^​ridge=(XTX+λI)−1XTy(3.44)

其中 I\mathbf{I}I 为 p×pp\times pp×p 的单位矩阵．注意到选择二次函数惩罚 βTβ\beta^T\betaβTβ，岭回归的解仍是 y\mathbf{y}y 的线性函数．解在求逆之前向矩阵 XTX\mathbf{X^TX}XTX 的对角元上加入正的常数值．即使 XTX\mathbf{X^TX}XTX 不是满秩，这样会使得问题非奇异，而且这是第一次将岭回归引入统计学中 (Hoerl and Kennard, 1970^{1）的主要动力．传统的岭回归的描述从定义 (3.44) 开始．我们选择通过 (3.41) 和 (3.42) 来阐述，因为这两式让我们看清楚了它是怎样实现的．}

图 3.8 展示了前列腺癌例子的岭回归系数估计，绘制成关于 df(λ)df(\lambda)df(λ) 的函数图象，df(λ)df(\lambda)df(λ) 为由惩罚 λ\lambdaλ 得到的有效自由度 (effective degrees of freedom)（由式 (3.50) 中定义）．在正交输入的情形下，岭回归估计仅仅是最小二乘估计的缩小版本，也就是 β^ridge=β^/(1+λ)\hat{\beta}^{ridge}=\hat{\beta}/(1+\lambda)β^​ridge=β^​/(1+λ)．

图 3.8 当惩罚参数 λ\lambdaλ 不同时，前列腺癌例子岭回归的变化曲线．画出系数关于有效自由度 df(λ)df(\lambda)df(λ) 的曲线．垂直直线画在 df=5.0df=5.0df=5.0 处，这是由交叉验证选择出来的．

当给定一个合适的先验分布，岭回归也可以从后验分布的均值或众数得到．具体地，假设 yi∼N(β0+xiTβ,σ2)y_i \sim N(\beta_0+x^T_i\beta,\sigma^2)yi​∼N(β0​+xiT​β,σ2)，参数 βj\beta_jβj​ 的分布均为 N(0,τ2)N(0,\tau^2)N(0,τ2)，每个都相互独立．则当 τ2\tau^2τ2 和 σ2\sigma^2σ2 值已知时，β\betaβ 后验分布密度函数的对数值（的负数）与 (3.41) 中花括号里面的表达式成比例**(weiya 注：原文直接说与花括号的表达式相等，但应该是常数倍)**，且 λ=σ2/τ2\lambda=\sigma^2/\tau^2λ=σ2/τ2（练习 3.6)．因此岭回归估计是后验分布的众数；又因分布为高斯分布，则也是后验分布的均值．

!!! info “weiya 注：Ex. 3.6”
将解答过程移至Issue 96: Ex. 3.6．

中心化输入矩阵 X\mathbf{X}X 的奇异值分解 (SVD) 让我们进一步了解了岭回归的本质．这个分解在许多统计方法分析中非常有用．N×pN\times pN×p 阶矩阵 X\mathbf{X}X 的 SVD 分解有如下形式
X=UDVT(3.45) \mathbf{X=UDV^T}\tag{3.45} X=UDVT(3.45)

这里 U\mathbf{U}U 和 V\mathbf{V}V 分别是 N×pN\times pN×p 和 p×pp\times pp×p 的正交矩阵，U\mathbf{U}U的列张成 XXX 的列空间，V\mathbf{V}V 的列张成 XXX 的行空间．D\mathbf{D}D 为 p×pp\times pp×p 的对角矩阵，对角元 d1≥d2≥⋯≥dp≥0d_1\ge d_2 \ge \cdots \ge d_p \ge 0d1​≥d2​≥⋯≥dp​≥0 称作 X\mathbf{X}X 的奇异值．如果一个或多个 dj=0d_j=0dj​=0，则 X\mathbf{X}X 为奇异的．

weiya 注：（以下关于奇异值分解内容来自张贤达的《矩阵分析与应用》）

奇异值分解最早由Beltrami在1873年对实正方矩阵提出来的．Beltrami从双线性函数
f(x,y)=xTAy,A∈Rn×m f(x,y)=x^TAy,\qquad A\in R^{n\times m} f(x,y)=xTAy,A∈Rn×m
出发，通过引入线性变换
x=Uξ,y=Vη x=U\xi,\qquad y=V\eta x=Uξ,y=Vη
将双线性函数变为
f(x,y)=ξTSη f(x,y)=\xi^TS\eta f(x,y)=ξTSη
式中
S=UTAV S=U^TAV S=UTAV
若选择UUU和VVV为正交矩阵，则他们的选择各存在n2−nn^2-nn2−n个自由度．他提出利用这些自由度使矩阵SSS的非对角元为0，即矩阵S=Σ=diag(σ1,σ2,…,σn)S=\Sigma=diag(\sigma_1,\sigma_2,\ldots,\sigma_n)S=Σ=diag(σ1​,σ2​,…,σn​)为对角矩阵．则
A=UΣVT A=U\Sigma V^T A=UΣVT
这是Beltrami于1873年得到的实正方矩阵的奇异值分解．后来，Autonne于1902年把奇异值分解推广到复正方矩阵；Eckart与Young于1939年又进一步把它推广到一般的长方形矩阵．因此，现在常将任意复长方矩阵奇异值分解定理称为Autonee-Eckart-Young定理，详见下述

令A∈Rm×nA\in R^{m\times n}A∈Rm×n(或Cm×nC^{m\times n}Cm×n),则存在正交（或酉）矩阵U∈Rm×mU\in R^{m\times m}U∈Rm×m(或Cm×mC^{m\times m}Cm×m)和V∈Rn×nV\in R^{n\times n}V∈Rn×n(或Cn×nC^{n\times n}Cn×n)使得
A=UΣVT(orUΣVH) A=U\Sigma V^T(or\quad U\Sigma V^H) A=UΣVT(orUΣVH)
式中
Σ=[Σ1OOO] \Sigma= \left[ \begin{array}{cc} \Sigma_1&O\\ O&O \end{array} \right] Σ=[Σ1​O​OO​]
且Σ1=diag(σ1,σ2,…,σr)\Sigma_1=diag(\sigma_1,\sigma_2,\ldots,\sigma_r)Σ1​=diag(σ1​,σ2​,…,σr​),其对角元素按照顺序
σ1>σ2⋯≥σr>0,r=rank(A) \sigma_1\gt \sigma_2\cdots\ge\sigma_r>0,\qquad r=rank(A) σ1​>σ2​⋯≥σr​>0,r=rank(A)
排列

以下说明来自wiki:
这里的奇异值分解称作Thin SVD
M=UnΣnVT M=U_n\Sigma_nV^T M=Un​Σn​VT
只有UUU中对应VTV^TVT中行向量的nnn个列向量被计算，UUU的其它列不被计算．UnU_nUn​为m×nm\times nm×n,Σn\Sigma_nΣn​为n×nn\times nn×n的对角矩阵，VVV是n×nn\times nn×n

利用奇异值分解，通过简化我们可以把最小二乘拟合向量写成

Xβ^ls=X(XTX)−1XTy=UUTy(3.46) \begin{aligned} \mathbf{X}\hat{\beta}^{ls}&=\mathbf{X(X^TX)^{-1}X^Ty}\\ &=\mathbf{UU^Ty}\tag{3.46} \end{aligned} Xβ^​ls​=X(XTX)−1XTy=UUTy​(3.46)

注意到 UTy\mathbf{U}^T\mathbf yUTy 是 y\mathbf{y}y 正交基 U\mathbf{U}U 下的坐标．同时注意其与 (3.33) 的相似性；Q\mathbf{Q}Q 和 U\mathbf{U}U 是 X\mathbf{X}X 列空间的两个不同的正交基（练习 3.8）．

!!! info “weiya 注：Ex. 3.8”
已解答，具体证明过程见 Issue 97: Ex. 3.8

现在岭回归的解为

Xβ^ridge=X(XTX+λI)−1XTy=UD(D2+λI)−1DUTy=∑j=1pujdj2dj2+λujTy(3.47) \begin{aligned} \mathbf{X}\hat{\beta}^{ridge}&=\mathbf{X}(\mathbf{X^TX}+\lambda \mathbf{I})^{-1}\mathbf{X^Ty}\\ &= \mathbf{UD}(\mathbf{D^2}+\lambda \mathbf{I})^{-1}\mathbf{DU^Ty}\\ &= \sum\limits_{j=1}^p\mathbf{u}_j\dfrac{d_j^2}{d_j^2+\lambda}\mathbf{u_j^Ty}\tag{3.47} \end{aligned} Xβ^​ridge​=X(XTX+λI)−1XTy=UD(D2+λI)−1DUTy=j=1∑p​uj​dj2​+λdj2​​ujT​y​(3.47)

其中 uj\mathbf{u}_juj​ 是 U\mathbf{U}U 的列向量．注意到因为 λ≥0\lambda \ge 0λ≥0，我们有 dj2/(dj2+λ)≤1d_j^2/(d^2_j+\lambda)\le 1dj2​/(dj2​+λ)≤1．类似线性回归，岭回归计算 y\mathbf{y}y 关于正规基 U\mathbf{U}U 的坐标．通过因子 dj2/(dj2+λ)d^2_j/(d^2_j+\lambda)dj2​/(dj2​+λ) 来收缩这些坐标．这意味着更小的 dj2d_j^2dj2​ 会在更大程度上收缩基向量的坐标．

dj2d_j^2dj2​ 值小意味着什么？中心化后的矩阵 X\mathbf{X}X 的奇异值分解是表示 X\mathbf{X}X 中主成分变量的另一种方式．样本协方差矩阵为 S=XTX/N\mathbf{S=X^TX}/NS=XTX/N

S=E((X−EX)T(X−EX))=XTX/N\mathbf{S={\color{red} E((X-EX)^T(X-EX))=}X^TX}/NS=E((X−EX)T(X−EX))=XTX/N，并且从 (3.45) 式我们得到
XTX=VD2VT(3.48) \mathbf{X^T X = VD^2V^T} \tag{3.48} XTX=VD2VT(3.48)

上式是 XTX\mathbf{X^TX}XTX（也是 SSS，取决于因子 NNN）的特征值分解 (eigen decomposition)．特征向量 vjv_jvj​（V\mathbf{V}V 的列向量）也称作 X\mathbf{X}X 的主成分（或 Karhunen-Loeve）方向．第一主成分方向 v1v_1v1​ 有下面性质：z1=Xv1\mathbf{z}_1=\mathbf{X}v_1z1​=Xv1​ 在所有 X\mathbf{X}X 列的标准化线性组合中有最大的样本方差．样本方差很容易看出来是

Var(z1)=Var(Xv1)=d12N(3.49) \rm{Var}(\mathbf{z}_1)=\rm{Var}(\mathbf{X}v_1)=\dfrac{d_1^2}{N}\tag{3.49} Var(z1​)=Var(Xv1​)=Nd12​​(3.49)

事实上 z1=Xv1=u1d1\mathbf{z}_1=\mathbf{X}v_1=\mathbf{u}_1d_1z1​=Xv1​=u1​d1​．导出变量 z1\mathbf{z_1}z1​ 称作 X\mathbf{X}X 的第一主成分，因此 u1\mathbf{u_1}u1​ 是标准化的第一主成分．后面的主成分 zjz_jzj​ 在与前一个保持正交的前提下有最大的方差 dj2/Nd_j^2/Ndj2​/N．所以，最后一个主成分有最小的方差．因此越小的奇异值 djd_jdj​ 对应 X\mathbf{X}X 列空间中方差越小的方向，并且岭回归在这些方向上收缩得最厉害．

图 3.9 展示了两个维度下部分数据点的主成分．如果我们考虑在这个区域（YYY 轴垂直纸面）内拟合线性曲面，数据的结构形态使得确定梯度时长方向会比短方向更精确．岭回归防止在短方向上估计梯度可能存在的高方差．隐含的假设是响应变量往往在高方差的输入方向上变化．这往往是个合理的假设，因为我们所研究的预测变量随响应变量变化而变化，而不需要保持不变．

图 3.9 部分输入数据点的主成分．最大主成分是使得投影数据方差最大的方向，最小主成分是使得方差最小的方向．岭回归将 y\mathbf{y}y 投射到这些成分上，然后对低方差成分的系数比高方差收缩得更厉害．

在图 3.7 中我们已经画了预测误差估计值关于 df(λ)df(\lambda)df(λ) 的曲线

df(λ)=tr[X(XTX+λI)−1XT]=tr(Hλ)=∑j=1pdj2dj2+λ(3.50) \begin{aligned} df(\lambda)&=tr[\mathbf{X}(\mathbf{X^TX}+\lambda\mathbf{I})^{-1}\mathbf{X}^T]\\ &=tr(\mathbf{H}_{\lambda})\\ &=\sum\limits_{j=1}^p\dfrac{d_j^2}{d_j^2+\lambda}\tag{3.50} \end{aligned} df(λ)​=tr[X(XTX+λI)−1XT]=tr(Hλ​)=j=1∑p​dj2​+λdj2​​​(3.50)

上面 λ\lambdaλ 的单调递减函数是岭回归拟合的有效自由度 (effective degrees of freedom)．通常在含 ppp 个变量的线性回归拟合中，拟合的自由度为 ppp，也就是无约束参数的个数．这里想法是尽管岭回归拟合中所有的 ppp 个系数都不为 0，但是它们在由 λ\lambdaλ 控制的约束下拟合．注意到当 λ=0\lambda=0λ=0（没有正则化）时 df(λ)=pdf(\lambda)=pdf(λ)=p，并且当 λ→∞\lambda\rightarrow \inftyλ→∞ 时 df(λ)→0df(\lambda)\rightarrow 0df(λ)→0．当然总是对于截距总有一个额外的自由度，事先 (apriori) 已经去掉了．这个定义将在 3.4.4 节和 7.4-7.6 节中详细介绍．图 3.7 中最小值在 df(λ)=5.0df(\lambda)=5.0df(λ)=5.0 处．表 3.3 表明岭回归将全最小二乘估计的测试误差降低了一小部分．

个人解读 拉格朗日乘子法

岭回归等式(3.41)其实就是对于(3.42)使用拉格朗日乘子法得到的，对于岭回归等式第一项其实几何意义上是一个椭圆形的抛物面，第二项几何意义上是一个圆柱形的抛物面，数学意义上第一项是最小二乘均方误差最小化，第二项是L2范数惩罚项。如果觉得前面的有关拉格朗日乘子法的介绍复杂，这篇博客中的简单介绍较易理解。

接下来原文给出了岭回归的参数解，其推导过程如下：
RSS(λ)=(y−Xβ)T(y−Xβ)+λβTβ=yTy−2βTXTy+βTXTXβ+λβTβ \begin{aligned} R S S(\lambda) &=(y-\mathbf{X} \beta)^{T}(y-\mathbf{X} \beta)+\lambda \beta^{T} \beta \\ &=y^{T} y-2 \beta^{T} \mathbf{X}^{T} y+\beta^{T} \mathbf{X}^{T} \mathbf{X} \beta+\lambda \beta^{T} \beta \end{aligned} RSS(λ)​=(y−Xβ)T(y−Xβ)+λβTβ=yTy−2βTXTy+βTXTXβ+λβTβ​
令误差RSS对β\betaβ的导数为0，则：
∇βRSS(λ)=−2XTy+2XTXβ+2λβ=0⇒(XTX+λI)β=XTy⇒β^ridge=(XTX+λI)−1XTy \begin{aligned} &\nabla_{\beta} R S S(\lambda)=-2 \mathbf{X}^{T} y+2 \mathbf{X}^{T} \mathbf{X} \beta+2 \lambda \beta=0\\ &\Rightarrow\left(\mathbf{X}^{T} \mathbf{X}+\lambda I\right) \beta=\mathbf{X}^{T} y\\ &\Rightarrow \hat{\beta}^{r i d g e}=\left(\mathbf{X}^{T} \mathbf{X}+\lambda I\right)^{-1} \mathbf{X}^{T} \mathbf{y} \end{aligned} ​∇β​RSS(λ)=−2XTy+2XTXβ+2λβ=0⇒(XTX+λI)β=XTy⇒β^​ridge=(XTX+λI)−1XTy​
由此得到(3.44)(3.44)(3.44)的结论。

作者在文中并没有定义什么是中心化输入矩阵XXX，中心化一个矩阵代表对于这个矩阵的列进行normalization，令其均值为0，方差为1，用的就是如下算法：

首先令
μj=1n∑i=1nxij \mu_{j}=\frac{1}{n} \sum_{i=1}^{n} x_{i j} μj​=n1​i=1∑n​xij​
然后对于每一个xj∈Xx_j \in Xxj​∈X，将其替换为xj−μjx_j-\mu_jxj​−μj​。之后令
σj2=1n−1∑i=1n(xij)2 \sigma_{j}^{2}=\frac{1}{n-1} \sum_{i=1}^{n}\left(x_{i j}\right)^{2} σj2​=n−11​i=1∑n​(xij​)2
再把每个xjx_jxj​替换为xjσj\frac{x_j}{\sigma_j}σj​xj​​。由此便可以对XXX进行奇异值分解了。

奇异值分解

许多人会有所疑问，一般的奇异值分解两边的矩阵都是方阵，中间的矩阵不是，就如下图所示：

这和书中所分解的结果是不一样的，书中的U不是方阵，中间的是方阵，命名这样的分解方法是Thin SVD，那么两者有什么区别吗？

实际上没有区别，假设U和Σ\SigmaΣ的乘积为B，矩阵B的每一列，对应的是矩阵U的线性组合的结果，对应的线性组合的系数，就是矩阵Σ\SigmaΣ的各列。但是，矩阵Σ\SigmaΣ的第N行以下全是0，所以线性组合的系数全是0，也就是U的第N列往后，对应的系数是0，他们根本没参与计算，所以可以直接丢掉，不影响结果。

那么再用奇异值分解之后的X=UDVTX=UDV^TX=UDVT（UUU是XXX的列空间，VVV是XXX的行空间）来重新表示岭回归的解的推导过程如下：
XTX=VDUTUDVT=VD2VT X^{T} X=V D U^{T} U D V^{T}=V D^{2} V^{T} XTX=VDUTUDVT=VD2VT
那么岭回归系数表达式则变化如下
β^ridge =(XTX+λI)−1XTy=(VD2VT+λVVT)−1VDUTy=(V(D2+λI)VT)−1VDUTy=V(D2+λI)−1DUTy \begin{aligned} \hat{\beta}^{\text {ridge }} &=\left(X^{T} X+\lambda I\right)^{-1} X^{T} \mathbf{y} \\ &=\left(V D^{2} V^{T}+\lambda V V^{T}\right)^{-1} V D U^{T} \mathbf{y} \\ &=\left(V\left(D^{2}+\lambda I\right) V^{T}\right)^{-1} V D U^{T} \mathbf{y} \\ &=V\left(D^{2}+\lambda I\right)^{-1} D U^{T} \mathbf{y} \end{aligned} β^​ridge ​=(XTX+λI)−1XTy=(VD2VT+λVVT)−1VDUTy=(V(D2+λI)VT)−1VDUTy=V(D2+λI)−1DUTy​
注意VVV是标准正交方阵，标准正交方阵的逆就是其转置，这就解释了上式倒数第二步如何变为最后一步。有关正交方阵的逆就是其转置证明如下：

因为是标准正交方阵，所以VTV=IV^TV=IVTV=I，这个等式证明很好理解，标准正交基的意思，就是空间中的单位正交基，各个基之间正交，转置乘以本身的结果中，对角元素就是各个基和本身相乘，也就是各个基向量的模，所以是1，非对角元素，就是各个基向量和别的基向量的乘积，相互正交，所以是0，所以结果是单位阵。那再在等式两边乘V−1V^{-1}V−1，就得最终结果
VT=V−1 V^T=V^{-1} VT=V−1
得到岭回归系数新的表达式之后，就可继续推得最终预测值为
y^ridge =Xβ^ridge =UD(D2+λI)−1DUTy \hat{y}^{\text {ridge }}=X \hat{\beta}^{\text {ridge }}=U D\left(D^{2}+\lambda I\right)^{-1} D U^{T} \mathbf{y} y^​ridge =Xβ^​ridge =UD(D2+λI)−1DUTy
注意D(D2+λI)−1DD\left(D^{2}+\lambda I\right)^{-1} DD(D2+λI)−1D是对角元素为dj2dj2+λ\dfrac{d_j^2}{d_j^2+\lambda}dj2​+λdj2​​的对角阵，令uju_juj​为UUU的列向量，可得最终结果
y^ridge=∑j=1pujdj2dj2+λujTy \hat y^{ridge}=\sum\limits_{j=1}^p\mathbf{u}_j\dfrac{d_j^2}{d_j^2+\lambda}\mathbf{u_j^Ty} y^​ridge=j=1∑p​uj​dj2​+λdj2​​ujT​y
此时看看最小二乘的对XXX进行奇异值分解代入后的预测值是什么形式
y^ls=Xβ^ls=UDVT(VD2VT)−1VDUTy=UDVT(V−TD−2V−1)VDUTy=UUTy=∑j=1puj(ujTy) \begin{aligned} \hat{y}^{\mathrm{ls}}=X \hat{\beta}^{\mathrm{ls}} &=U D V^{T}\left(V D^{2} V^{T}\right)^{-1} V D U^{T} \mathrm{y} \\ &=U D V^{T}\left(V^{-T} D^{-2} V^{-1}\right) V D U^{T} \mathrm{y} \\ &=U U^{T} \mathrm{y} \\ &=\sum_{j=1}^{p} u_{j}\left(u_{j}^{T} \mathrm{y}\right) \end{aligned} y^​ls=Xβ^​ls​=UDVT(VD2VT)−1VDUTy=UDVT(V−TD−2V−1)VDUTy=UUTy=j=1∑p​uj​(ujT​y)​
由此可见，岭回归中，ujTyu_{j}^{T} \mathrm{y}ujT​y的内积现在被dj2dj2+λ\dfrac{d_j^2}{d_j^2+\lambda}dj2​+λdj2​​收缩了。

所以这就印证了越大的λ\lambdaλ或越小的djd_jdj​会更大程度上收缩基向量的坐标。

之后就提及了特征值分解，有关特征值分解和奇异值分解的几何意义可以参考这篇博客，更加有助于直观理解，但博客中的内容和书中的理解还仍有差异，而这篇博客中提及的概念同样有助浅层理解，想形成属于自己体系的理解还是跟着本书的思想走。

特征值分解过后，对于第一主成分方向上样本方差的表达式证明如下
Var⁡(Xvi)=viTVar⁡(X)vi=viT(XTXN)vi=viTVD2VTvi/N \begin{aligned} \operatorname{Var}\left(\mathbf{X} v_{i}\right)=& v_{i}^{T} \operatorname{Var}(\mathbf{X}) v_{i} \\ =& v_{i}^{T}\left(\frac{\mathbf{X}^{T} \mathbf{X}}{N}\right) v_{i} \\ =& v_{i}^{T} \mathbf{V D}^{2} \mathbf{V}^{T} v_{i} / N \end{aligned} Var(Xvi​)===​viT​Var(X)vi​viT​(NXTX​)vi​viT​VD2VTvi​/N​
由于V是标准正交矩阵，非对角元素相乘，就是各个基向量和别的基向量的乘积，相互正交，所以是0。因此上式就
=viTvidi2vitvi/N =v_{i}^{T} v_{i} d_{i}^{2} v_{i}^{t} v_{i} / N =viT​vi​di2​vit​vi​/N
对角元素相乘就是各个基和本身相乘，也就是各个基向量的模，所以是1，所以最终可得
Var⁡(Xvi)=di2N \operatorname{Var}\left(\mathbf{X} v_{i}\right)=\frac{d_{i}^{2}}{N} Var(Xvi​)=Ndi2​​
由此引出重要结论：越小的奇异值 djd_jdj​ 对应 X\mathbf{X}X 列空间中方差越小的方向，并且岭回归在这些方向上收缩得最厉害。

接下来就要讲到Lasso回归了，岭回归与Lasso回归最终都实现了约束参数防止过拟合的效果，但Lasso可将作用小的特征的参数直接训练成0，从而实现了特征筛选降维的操作，因此Lasso回归也被称为是温和的连续子集选择。

Lasso

!!! note “weiya 注：”
lasso 是 “Least absolute shrinkage and seleetion operator” 的首字母缩写．

lasso 像岭回归一样是个收缩方法，有微妙但很重要的区别．lasso 估计定义如下

β^lasso=arg⁡min⁡β∑i=1N(yi−β0−∑j=1pxijβj)2subject to ∑j=1p∣βj∣≤t(3.51) \begin{aligned} \hat{\beta}^{lasso}&=\underset{\beta}{\arg\min}\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2\\ &\text{subject to }\sum\limits_{j=1}^p\vert\beta_j\vert\le t \tag{3.51} \end{aligned} β^​lasso​=βargmin​i=1∑N​(yi​−β0​−j=1∑p​xij​βj​)2subject to j=1∑p​∣βj​∣≤t​(3.51)

正如在岭回归中一样，我们可以通过标准化预测变量来对常数 β0\beta_0β0​ 再参量化；β^0\hat{\beta}_0β^​0​ 的解为 yˉ\bar{y}yˉ​，并且后面我们拟合无截距的模型（练习 3.5）．

!!! info “weiya 注：Ex. 3.5”
已解答，详细证明过程见 Issue 95: Ex. 3.5

在信号处理中，lasso 也被称作 basis pursuit (Chen et al., 1998²⁾

我们也可以把 lasso 问题等价地写成拉格朗日形式 (Lagrangian form)

β^lasso=arg⁡min⁡β{∑i=1N(yi−β0−∑j=1pxijβj)2+λ∑j=1p∣βj∣}(3.52) \hat{\beta}^{lasso}=\underset{\beta}{\arg\min}\Big\{\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2+\lambda\sum\limits_{j=1}^p\vert\beta_j\vert\Big\}\tag{3.52} β^​lasso=βargmin​{i=1∑N​(yi​−β0​−j=1∑p​xij​βj​)2+λj=1∑p​∣βj​∣}(3.52)

注意到这与岭回归问题 (3.42) 或 (3.41) 的相似性：L2L_2L2​ 的岭回归惩罚 ∑1pβj2\sum_1^p\beta^2_j∑1p​βj2​ 替换为 L1L_1L1​ 的 lasso 惩罚∑1p∣βj∣\sum_1^p\vert\beta_j\vert∑1p​∣βj​∣．后一约束使得解在 yiy_iyi​ 处非线性，并且在岭回归中没有相近的表达式．计算 lasso 的解是一个二次规划问题，尽管我们在 3.4.4 节看到当 λ\lambdaλ 不同时计算解的整个路径存在与岭回归同样计算量的有效算法．由于该约束的本质，令 ttt 充分小会造成一些参数恰恰等于 0．因此 lasso 完成一个温和的连续子集选择．如果所选的 ttt 大于t0=∑1p∣β^j∣t_0=\sum_1^p\vert\hat{\beta}_j\vertt0​=∑1p​∣β^​j​∣（其中 β^j=β^jls\hat{\beta}_j=\hat{\beta}_j^{ls}β^​j​=β^​jls​，β^jls\hat{\beta}_j^{ls}β^​jls​ 为最小二乘估计），则 lasso 估计为 β^j\hat{\beta}_jβ^​j​．另一方面，当 t=t0/2t=t_0/2t=t0​/2，最小二乘系数平均收缩 50%50\%50%．然而，收缩的本质不是很显然，我们将在 3.4.4 节进一步研究．类似在变量子集选择中子集的大小，或者岭回归的惩罚参数，应该自适应地选择 ttt 使预测误差期望值的估计最小化．

图 3.7 中，为了方便解释，我们已经画出 lasso 的预测误差估计关于标准化参数 s=t/∑1p∣β^j∣s=t/\sum^p_1\vert\hat{\beta}_j\verts=t/∑1p​∣β^​j​∣ 的曲线．通过 10 折交叉验证选择 s≈0.36s\approx 0.36s≈0.36；这使得 4 个系数为 0（表 3.3 的第 5 列）．最终模型有第二低的测试误差，比全最小二乘模型略低，但是测试误差估计的标准误差（表 3.3 的最后一行）相当大．

图 3.10 显示了当惩罚参数 s=t/∑1p∣β^j∣s=t/\sum_1^p\vert\hat{\beta}_j\verts=t/∑1p​∣β^​j​∣ 不同时的 lasso 系数．当 s=1.0s=1.0s=1.0 时为最小二乘估计；当 s→0s\rightarrow 0s→0 时下降为 0．该下降不总是严格单调的，尽管例子中确实是．在 s=0.36s=0.36s=0.36 处画了垂直直线，该值通过交叉验证来选择．

图 3.10 当惩罚参数 ttt 变化时的 lasso 系数曲线．图中画了系数关于 s=t/∑1p∣β^j∣s=t/\sum^p_1\vert\hat{\beta}_j\verts=t/∑1p​∣β^​j​∣ 的曲线．垂直直线画在 s=0.36s=0.36s=0.36 处，该值通过交叉验证来选择．比较 65 页的图 3.8，lasso 曲线会达到 0，然而岭回归不会．曲线是分段线性的，所以只计算显示点处的值；详见 3.4.4 节．

讨论：子集的选择，岭回归，Lasso

这部分我们讨论并且比较至今为止有约束的线性回归模型的三种方法：子集选择、岭回归和 lasso．

在正交输入矩阵的情况下，三种过程都有显式解．每种方法对最小二乘估计 β^j\hat{\beta}_jβ^​j​ 应用简单的变换，详见表 3.4．

表 3.4 在 X\mathbf{X}X 为正规列情形下 βj\beta_jβj​ 的估计值．MMM 和 λ\lambdaλ 是通过对应的手段选择的常数；符号标记变量的符号（±1\pm 1±1），而且 x+x_+x+​ 记 xxx 的正数部分．下面的表格中，估计值由红色虚线来显示．灰色的 45∘45^{\circ}45∘ 直线作为参照显示了无约束的估计．

岭回归做等比例的收缩．lasso 通过常数因子 λ\lambdaλ 变换每个系数，在 0 处截去．这也称作“软阈限”，而且用在 5.9 节中基于小波光滑的内容中．最优子集选择删掉所有系数小于第 MMM 个大系数的变量；这是“硬阈限”的一种形式．

回到非正交的情形，一些图象可以帮助了解它们之间的关系．当只有两个参数时图 3.11 描绘了 lasso（左）和岭回归（右）．残差平方和为椭圆形的等高线，以全最小二乘估计为中心．岭回归的约束区域为圆盘 β12+β22≤t\beta_1^2+\beta_2^2\le tβ12​+β22​≤t，lasso 的约束区域为菱形∣β1∣+∣β2∣≤t\vert\beta_1\vert+\vert\beta_2\vert\le t∣β1​∣+∣β2​∣≤t．两种方式都寻找当椭圆等高线达到约束区域的第一个点．与圆盘不同，菱形 (diamond) 有角；如果解出现在角上，则有一个参数 βj\beta_jβj​ 等于 0．当 p>2p > 2p>2，菱形变成了偏菱形 (rhomboid)，而且有许多角，平坦的边和面；对于参数估计有更多的可能为 0．

图 3.11 lasso (左)和岭回归（右）的估计图象．图中显示了误差的等高线和约束函数．实心蓝色区域分别为约束区域∣β1∣+∣β2∣≤t\vert\beta_1\vert+\vert\beta_2\vert\le t∣β1​∣+∣β2​∣≤t以及β12+β22≤t2\beta^2_1+\beta_2^2\le t^2β12​+β22​≤t2，红色椭圆为最小二乘误差函数的等高线．

我们可以把岭回归和 lasso 一般化，并且可以看成是贝叶斯估计．考虑下面准则
β~=arg⁡min⁡β{∑i=1N(yi−β0−∑j=1pxijβj)2+λ∑j=1p∣βj∣q}(3.53) \tilde{\beta}=\underset{\beta}{\arg\min}\Big\{\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2+\lambda\sum\limits_{j=1}^p\vert\beta_j\vert^q\Big\}\tag{3.53} β~​=βargmin​{i=1∑N​(yi​−β0​−j=1∑p​xij​βj​)2+λj=1∑p​∣βj​∣q}(3.53)

其中 q≥0q\ge 0q≥0．图 3.12 显示了两个输入情形下常数值 ∑j∣βj∣q\sum_j\vert\beta_j\vert^q∑j​∣βj​∣q 的等高线．

图 3.12 给定值 qqq 下常数值 ∑j∣βj∣q\sum_j\vert\beta_j\vert^q∑j​∣βj​∣q 的等高线．

将 ∣βj∣q\vert\beta_j\vert^q∣βj​∣q 看成 βj\beta_jβj​ 的先验概率密度的对数值，同样有参数先验分布的等高线．q=0q=0q=0 对应变量子集选择，惩罚项是简单地统计非零参数的个数；q=1q=1q=1 对应 lasso，q=2q=2q=2 对应岭回归．注意到 q≤1q\le 1q≤1，先验在各方向上不是均匀的，而是更多地集中在坐标方向上．对应 q=1q=1q=1 情形的先验分布是关于每个输入变量是的独立的二重指数分布（或者 Laplace 分布），概率密度为(1/2τ)exp(−∣β∣)/τ(1/2\tau)exp(-\vert\beta\vert)/\tau(1/2τ)exp(−∣β∣)/τ 并且 τ=1/λ\tau=1/\lambdaτ=1/λ．q=1q=1q=1 的情形（lasso）是使得约束区域为凸的最小 qqq 值；非凸约束区域使得优化问题很困难．

从这点看，lasso、岭回归和最优子集选择是有着不同先验分布的贝叶斯估计．然而，注意到它们取自后验分布的众数，即最大化后验分布．在贝叶斯估计中使用后验分布的均值更加常见．岭回归同样是后验分布的均值，但是 lasso 和最优子集选择不是．

再一次观察准则 (3.53)，我们可能尝试除 0，1，2 外的其它 qqq 值．尽管有人可能从数据中估计 qqq，我们的经验表明引入额外的方差不值得．q∈(1,2)q\in (1,2)q∈(1,2) 表明在 lasso 和岭回归之间进行权衡．当 q>1q > 1q>1 时尽管 ∣βj∣q\vert\beta_j\vert^q∣βj​∣q 在 0 处可导，但是并没有lasso（q=1q=1q=1）的令系数恰巧为零的性质．部分由于这个原因并且考虑计算易处理，Zou and Hastie (2005)^{3 引入弹性惩罚}

λ∑j=1p(αβj2+(1−α)∣βj∣)(3.54) \lambda \sum\limits_{j=1}^p(\alpha\beta_j^2+(1-\alpha)\vert\beta_j\vert)\tag{3.54} λj=1∑p​(αβj2​+(1−α)∣βj​∣)(3.54)

这是一种岭回归和 lasso之间的平衡．图 3.13 比较了 q=1.2q=1.2q=1.2 下的 LqL_qLq​ 惩罚以及 α=0.2\alpha=0.2α=0.2 的弹性网惩罚；很难从肉眼来观察出差异．弹性网像 lasso 一样选择变量，同时像岭回归一样收缩相关变量的系数．同时考虑了 LqL_qLq​ 惩罚的计算优势．我们将在 18.4节介绍弹性网惩罚．

图3.13 q=1.2q=1.2q=1.2 时 ∑j∣βj∣q\sum_j\vert\beta_j\vert^q∑j​∣βj​∣q 为常数值的轮廓线（左图）以及 α=0.2\alpha=0.2α=0.2 时弹性网惩罚 ∑j(αβj2+(1−α)∣βj∣)\sum_j(\alpha\beta_j^2+(1-\alpha)\vert\beta_j\vert)∑j​(αβj2​+(1−α)∣βj​∣) 为常数值的轮廓线（右图）．尽管看起来很相似，弹性网有尖角（不可导），而 q=1.2q=1.2q=1.2 的惩罚不会有尖角．

最小角回归

最小角回归 (LAR) 是相对较新的方法 (Efron et al., 2004^{4)，而且可以看成是一种向前逐步回归（3.3.2 节）的“民主 (democratic)”版本．正如我们将看到的，LAR 与 lasso 联系紧密，并且事实上提供了如图 3.10 所示的计算整个 lasso 路径的非常有效的算法．}

!!! note “weiya 注：”
在 Efron 的论文中，最小角回归缩写为 LARS，我们这里仍以 ESL 书上的缩写为准——LAR

向前逐步回归逐步建立模型，每次添加一个变量．每一步，它选出最好的变量加入活跃集，然后更新最小二乘来加入所有的活跃变量．

最小角回归采用类似的策略，但是仅仅加入一个变量应有的程度．第一步它确定与响应变量最相关的变量．不是完全的拟合该变量，LAR 使得该变量的系数向最小二乘值连续变化（使得它与进化的残差之间的相关系数绝对值降低）．只要其他变量与残差的相关性与该变量和残差的相关性相等，则该过程暂停．第二个变量加入活跃集，然后它们的系数一起以保持相关性相等并降低的方式变化．这个过程一直继续直到所有的变量都在模型中，然后在全最小二乘拟合处停止．算法 3.2 给出了详细过程．第 5 步的终止条件需要一些解释．如果 p>N−1p >N-1p>N−1，LAR 算法经过 N−1N-1N−1 步达到 0 残差解（−1-1−1 是因为我们已经对数据进行了中心化）

算法 3.2 最小角回归

对预测变量进行标准化处理得到零均值和单位范数．以残差向量 r=y−yˉ,β1,…,βp=0\mathbf{r=y-\bar{y}},\beta_1,\ldots,\beta_p=0r=y−yˉ​,β1​,…,βp​=0 开始． 找出与 r\mathbf{r}r 最相关的预测变量 xj\mathbf x_jxj​ 从 0 开始移动 βj\beta_jβj​ 一直到最小二乘系数 ⟨xj,r⟩\langle\mathbf x_j, \mathbf r\rangle⟨xj​,r⟩，直到存在其它的预测变量 xk\mathbf x_kxk​ 使得其与当前残差的相关性等于 xj\mathbf x_jxj​ 与当前残差的相关性．（Hytn注：相关性相等几何意义上就是角平分线） 再由当前残差在 (xj,xk)(\mathbf x_j,\mathbf x_k)(xj​,xk​) 上的联合最小二乘系数方向上移动 βj\beta_jβj​ 和 βk\beta_kβk​，直到存在其它的预测变量 xlx_lxl​ 与当前残差的相关性和当前残差与 (xj,xk)(\mathbf x_j,\mathbf x_k)(xj​,xk​) 的相关性相等． 按这种方式继续直到所有的 ppp 个预测变量加入到模型中．经过 min⁡(N−1,p)\min(N-1, p)min(N−1,p) 步，我们达到了全最小二乘的解．

!!! tip “weiya 注：LAR 示意图”
结合原论文的示意图能帮助理解最小角回归的逻辑．

(Hytn注：就是不断尝试以最快的方向迅速拟合到预测变量，图中预测变量和x1x_1x1​的夹角小，所以在x1x_1x1​方向上前进一段距离，令当前点和预测变量间的向量μ2\mu_2μ2​形成二等分线，)

假设 Ak\mathcal A_kAk​ 是第 kkk 步开始时的变量活跃集，βAk\beta_{\mathcal A_k}βAk​​ 是这一步中变量的系数向量；其中有 k−1k-1k−1 个非零值，刚刚进入的变量系数值为 0．如果当前残差为 rk=y−XAkβAk\mathbf r_k=\mathbf y-\mathbf X_{\mathcal A_k}\beta_{\mathcal A_k}rk​=y−XAk​​βAk​​，则当前步的方向为
δk=(XAkTXAk)−1XAkTrk(3.55) \delta_k=(\mathbf X^T_{ \mathcal A_k}\mathbf X_{\mathcal A_k})^{-1}\mathbf X^T_{\mathcal A_k}\mathbf r_k \tag{3.55} δk​=(XAk​T​XAk​​)−1XAk​T​rk​(3.55)

然后系数迭代为 βAk(α)=βAk+α⋅δk\beta_{\mathcal A_k} (\alpha) = \beta_{\mathcal A_k} + \alpha · \delta_kβAk​​(α)=βAk​​+α⋅δk​．练习 3.23 证明这种方式下选择的方向满足断言：保持（各个预测变量与残差间的）相关系数相等和递减（tied and decreasing）．

!!! info “weiya 注：Ex. 3.23”
已解决，具体证明过程参见Issue 100: Ex. 3.23．起初翻译时，对 tied 的理解不够，通过求解该练习题，认为 tied 意思其实就是各个预测变量与残差之间的相关系数保持相等．

如果该步的开始拟合向量为 f^k\hat{\mathbf f}_kf^k​，则迭代为 f^k(α)=fk+α⋅uk\hat{\mathbf f}_k(\alpha)=\mathbf f_k+\alpha\cdot\mathbf u_kf^k​(α)=fk​+α⋅uk​，其中 uk=XAkδk\mathbf u_k=\mathbf X_{\mathcal A_k}\delta_kuk​=XAk​​δk​ 是新的拟合方向．“最小角”由该过程的几何解释得到；uk\mathbf u_kuk​ 使得活跃集 Ak{\mathcal A}_kAk​中预测变量间的角度最小（练习 3.24）．

!!! info “weiya 注：Ex. 3.24”
已解决，详见 Issue 101: Ex. 3.24，欢迎交流讨论．

图 3.14 使用模拟数据显示了相关系数的绝对值下降以及每一步 LAR 算法中变量进入的顺序．

图 3.14：通过 6 个预测变量的拟合数据集，每一步 LAR 过程中的相关性绝对值的变化．图象上方的标签表示在每一步哪些变量加进了活跃集．步长是用单位 L1L_1L1​ 弧长来测量的．

由构造知 LAR 的系数以一种分段线性的方式进行改变．图 3.15（左图）显示了 LAR 系数曲线作为 L1L_1L1​ 弧长的函数曲线．

!!! note “weiya 注：原书脚注”
L1L_1L1​ arc length：可导曲线 β(s),s∈[0,S]\beta(s), s \in [0,S]β(s),s∈[0,S] 的 L1L_1L1​ 弧长为 TV(β,S)=∫0S∥β˙(s)∥1dsTV(\beta,S)=\int_0^S\Vert\dot{\beta}(s)\Vert_1dsTV(β,S)=∫0S​∥β˙​(s)∥1​ds，其中 β˙(s)=∂β(s)/∂s\dot{\beta}(s)=\partial\beta(s)/\partial sβ˙​(s)=∂β(s)/∂s．对于分段 LAR 函数曲线，这相当于从这一步到下一步系数的 L1L_1L1​ 范数变化之和．

图 3.15：左图显示了 LAR 系数作为 L1L_1L1​ 长度的函数在模拟数据上的图象．右图显示了 Lasso 的图象．它们大概在 L1L_1L1​ 弧长为 18 之前（深蓝色的系数曲线通过 0）都是完全相同的.

注意到我们不需要走很小的步以及重新检查步骤 3 的相关系数；应用预测变量的协方差和算法的分段线性性质，我们可以在每一步开始计算出确切的步长（练习 3.25）．

!!! info “weiya 注：Ex. 3.25”
已解决，详见 Issue 98: Ex. 3.25，欢迎讨论交流！

图 3.15 的右图展示了对同样数据的 lasso 系数曲线．几乎与左图相同，当绿色曲线通过 0 时首次出现不同．对于前列腺癌数据，LAR 系数曲线显示与图 3.10 的 lasso 曲线相同，该曲线从不经过 0．这些观测值促使对 LAR 算法进行简单修改，给出了整个 lasso 路径，它同样也是分段线性的．

算法 3.2a 最小角回归：Lasso修正

4a. 如果一个非零的系数达到0，则从变量的活跃集中删除该变量并且重新计算当前的联合最小二乘方向．

LAR(lasso) 算法是非常有效的，需要用 ppp 个预测变量的单最小二乘拟合的相同步骤进行计算．最小角回归总是需要 ppp 步达到全最小二乘估计．lasso 路径可能超过 ppp 步，尽管这两者经常是非常相似的．经过 lasso 修正的 3.2a 的算法 3.2 是计算任何一个lasso 问题的有效方式，特别是当 p>>Np > >Np>>N．Osborne et al. (2000a)^{5 也发现了计算 lasso 的分段线性的路径，他们称之为同伦 (homotopy) 算法．}

我们已经给出一个为什么这些过程很相似的启发式的论据．尽管 LAR 算法是用相关性来叙述的，但如果输入特征是标准化的，它与内积是等价的并且用内积更简单．假设 A\mathcal AA 是算法中某些步的变量活跃集，它们与当前残差 y−Xβ\mathbf y -\mathbf X\betay−Xβ 的内积的绝对值是结合在一起的．我们可以表达成

xjT(y−Xβ)=γ⋅sj,∀j∈A(3.56) \mathbf x_j^T(\mathbf y-\mathbf X\beta)=\gamma\cdot s_j,\forall j\in {\mathcal A} \tag{3.56} xjT​(y−Xβ)=γ⋅sj​,∀j∈A(3.56)

其中 sj∈{−1,1}s_j\in\{-1,1\}sj​∈{−1,1} 表示内积的符号，γ\gammaγ 是普通的数值．并且 ∣xkT(y−Xβ)∣≤γ  ∀k∉A\vert \mathbf x_k^T(\mathbf y-\mathbf X\beta)\vert\le \gamma\; \forall k\notin \mathcal A∣xkT​(y−Xβ)∣≤γ∀k∈/​A．现在我们考虑 (3.52) 的 lasso 准则，我们可以写成向量形式

R(β)=12∥y−Xβ∥22+λ∥β∥1(3.57) R(\beta)=\frac{1}{2}\Vert\mathbf y-\mathbf X\beta\Vert_2^2+\lambda\Vert\beta\Vert_1\tag{3.57} R(β)=21​∥y−Xβ∥22​+λ∥β∥1​(3.57)

令 B\mathcal BB 为在给定 λ\lambdaλ 值下解中的变量的活跃集．对于这些变量 R(β)R(\beta)R(β) 是可导的，并且平稳条件 (stationary condition) 为(Hytn注：也就是活跃变量导数为0的时候)

xjT(y−Xβ)=λ⋅sign(βj),∀j∈B(3.58) \mathbf x_j^T(\mathbf y-\mathbf X\beta)=\lambda\cdot \rm{sign}(\beta_j),\forall j\in {\mathcal B}\tag{3.58} xjT​(y−Xβ)=λ⋅sign(βj​),∀j∈B(3.58)

比较 3.58 和 3.56 ，我们看到只有当 βj\beta_jβj​ 的符号与内积的符号相同时，这两个等式才相同．这也就是为什么 LAR 算法和 lasso 当一个活跃系数经过零开始出现不同；对于不满足条件 3.58 的变量，会被踢出活跃集 B\mathcal BB．练习 3.23 证明了这些等式表明随λ\lambdaλ 减小的分段线性系数曲线．对于不活跃的变量的平稳条件要求

∣xkT(y−Xβ)∣≤λ,∀k∉B(3.59) \vert\mathbf x_k^T(\mathbf y-\mathbf X\beta)\vert\le\lambda,\forall k\notin{\mathcal B}\tag{3.59} ∣xkT​(y−Xβ)∣≤λ,∀k∈/​B(3.59)

这与 LAR 算法一致．

图 3.16：LAR、lasso、向前逐步、向前逐渐（FS）和增长向前逐渐（FS0FS_0FS0​）回归之间的比较．设定与图3.6相同，除了这里N=100N=100N=100而不是300.这里较慢的FSFSFS回归最终表现得比向前逐步好．LAR和lasso表现得和FS、FS0FS_0FS0​相似．因为这些过程采取不同的步数（根据模拟复制和方法），我们画出最小二乘拟合的MSE关于整体L1L_1L1​弧长的片段的函数．

图 3.16 将 LAR 和 lasso 与向前逐步（forward stepwise）和向前逐渐（forward stagewise）回归进行比较．设定与图 3.6 是相同的，除了这里的 N=100N=100N=100 而不是 300300300，所以这个问题更加困难．我们可以看到增长性更快的向前逐步很快地过拟合（10 个变量加入模型中之前是很好的），最终比增长性较慢的向前逐渐回归表现得更差．LAR 和 lasso 的行为与向前逐渐回归相似．增长的向前逐渐回归与 LAR 和 lasso 类似，并且将在 3.8.1 节中描述．

个人解读

书中对于最小角回归的描述如下

首先对于所有的predictors做一个中心化处理（均值0，单位方差），开始的时候所有的回归系数都是0，除了截距以外，刚开始的残差就是r=y−yˉr=y-\bar yr=y−yˉ​，β0=yˉ\beta_0=\bar yβ0​=yˉ​。

现在k=1进行第一步，找到和当前残差r1=y−yˉr_1=y-\bar yr1​=y−yˉ​最相关的 predictor xjx_jxj​。此时就有激活的方向A1={xj}\mathcal{A}_{1}=\left\{x_{j}\right\}A1​={xj​}和激活的系数βA1=[0]\beta_{\mathcal{A}_{1}}=[0]βA1​​=[0]。

把βj\beta_jβj​从初值0向着该值移动
δ1=(XA1TXA1)−1XA1Tr1=xjTr1xjTxj=xjTr1 \delta_{1}=\left(X_{\mathcal{A}_{1}}^{T} X_{\mathcal{A}_{1}}\right)^{-1} X_{\mathcal{A}_{1}}^{T} r_{1}=\frac{x_{j}^{T} r_{1}}{x_{j}^{T} x_{j}}=x_{j}^{T} r_{1} δ1​=(XA1​T​XA1​​)−1XA1​T​r1​=xjT​xj​xjT​r1​​=xjT​r1​
注意这个值的几何意义，就是残差r1r_1r1​在最相关方向xjx_jxj​上面的投影，也就是回归系数的最小二乘值，为什么不直接移动满就是算法的特性，停止条件就是出现另一个和残差的相关性跟当前predictor一样的predictor。那么移动到底多少，就有如下表达式说明
βA1(α)≡βA1+αδ1=0+αxjTr1=(xjTr1)α for 0≤α≤1 \beta_{\mathcal{A}_{1}}(\alpha) \equiv \beta_{\mathcal{A}_{1}}+\alpha \delta_{1}=0+\alpha x_{j}^{T} r_{1}=\left(x_{j}^{T} r_{1}\right) \alpha \quad \text { for } \quad 0 \leq \alpha \leq 1 βA1​​(α)≡βA1​​+αδ1​=0+αxjT​r1​=(xjT​r1​)α for 0≤α≤1
得出回归系数之后实际移动的路径（也代表已经拟合好的值）就为
f^1(α)=XA1βA1(α)=(xjTr1)αxj \hat{f}_{1}(\alpha)=X_{\mathcal{A}_{1}} \beta_{\mathcal{A}_{1}}(\alpha)=\left(x_{j}^{T} r_{1}\right) \alpha x_{j} f^​1​(α)=XA1​​βA1​​(α)=(xjT​r1​)αxj​
那么还剩下的残差就是
r(α)=y−yˉ−α(xjTr1)xj=r1−α(xjTr1)xj r(\alpha)=y-\bar{y}-\alpha\left(x_{j}^{T} r_{1}\right) x_{j}=r_{1}-\alpha\left(x_{j}^{T} r_{1}\right) x_{j} r(α)=y−yˉ​−α(xjT​r1​)xj​=r1​−α(xjT​r1​)xj​
那么当前predictor和这个残差的联系就可描述为
xjT(r1−α(xjTr1)xj)=xjTr1−α(xjTr1)=(1−α)xjTr1 x_{j}^{T}\left(r_{1}-\alpha\left(x_{j}^{T} r_{1}\right) x_{j}\right)=x_{j}^{T} r_{1}-\alpha\left(x_{j}^{T} r_{1}\right)=(1-\alpha) x_{j}^{T} r_{1} xjT​(r1​−α(xjT​r1​)xj​)=xjT​r1​−α(xjT​r1​)=(1−α)xjT​r1​
而其他predictor例如xkx_kxk​和这个残差的联系就是
xkT(r1−α(xjTr1)xj)=xkTr1−α(xjTr1)xkTxj x_{k}^{T}\left(r_{1}-\alpha\left(x_{j}^{T} r_{1}\right) x_{j}\right)=x_{k}^{T} r_{1}-\alpha\left(x_{j}^{T} r_{1}\right) x_{k}^{T} x_{j} xkT​(r1​−α(xjT​r1​)xj​)=xkT​r1​−α(xjT​r1​)xkT​xj​

LAR 和 Lasso 自由度公式

假设我们通过最小角回归过程拟合了线性模型，在某步 k<pk < pk

ttt 得到约束情况下的全最小二乘拟合．我们需要多少参数，或者自由度？

首先考虑采用 kkk 个特征的子集的线性回归．如果这个子集是没有通过训练数据而事先确定好的，然后在该拟合模型中的自由度定义为kkk．当然，在经典统计学中，线性独立参数的个数也就是自由度．另外地，假设我们用一个最优子集选择确定了最优的 kkk 个预测变量．于是结果模型中有 kkk 个参数，但在某种意义上我们用了大于 kkk 个的自由度．

我们需要一个对于自适应拟合模型的有效自由度的一般定义．我们定义拟合向量 y^=(y^1,y^2,…,y^N)\hat{\mathbf y}=(\hat y_1,\hat y_2,\ldots,\hat y_N)y^​=(y^​1​,y^​2​,…,y^​N​) 的自由度为
df(y^)=1σ2∑i=1NCov(y^i,yi)(3.60) \rm{df}(\hat{\mathbf y})=\frac{1}{\sigma^2}\sum\limits_{i=1}^N\rm{Cov}(\hat y_i,y_i)\tag{3.60} df(y^​)=σ21​i=1∑N​Cov(y^​i​,yi​)(3.60)
这里 Cov(y^i,yi)\rm{Cov}(\hat y_i,y_i)Cov(y^​i​,yi​) 指的是预测值 y^i\hat y_iy^​i​ 和其对应的输出向量 yiy_iyi​ 之间的协方差．直观上看有意义：当拟合数据越困难，协方差会越大，从而 df(y^)\rm{df}(\hat{\mathbf y})df(y^​) 越大．表达式（3.60）是一个有用的自由度的概念，可以应用到任何模型的预测向量 y^\hat{\mathbf y}y^​．其中包括那些对训练数据自适应拟合的模型．这个定义将在 7.4-7.6 节 中进一步讨论．

现在对于有 kkk 个固定预测变量的线性回归模型，可以简单地证明 df(y^)=k\rm{df}(\hat{\mathbf y})=kdf(y^​)=k．同样地，对于岭回归，这一定义导出表达式（3.50）的闭型解 (closed-form)：df(y^)=tr(Sλ)\rm{df}(\hat{\mathbf{y}})=\rm{tr}(\mathbf S_\lambda)df(y^​)=tr(Sλ​)．

!!! note “weiya 注：closed form expression”
根据 wiki: Closed-form expression，closed-form 是指可以进行有限次赋值的表达式，依此理解，有显示解的为 closed form．

在这些情况下，(3.60) 可以很简单地进行赋值因为 y^=Hλy\hat{\mathbf{y}}=\mathbf{H}_\lambda\mathbf yy^​=Hλ​y 关于 y\mathbf yy 是线性的．如果我们考虑在大小为 kkk 的最优子集选择中的定义 (3.60)，似乎显然有 df(y^)\rm{df}(\hat{\mathbf y})df(y^​) 会大于 kkk，并且可以通过运用模拟的方法直接地估计 Cov(y^i,yi)/σ2\rm{Cov}(\hat y_i,y_i)/\sigma^2Cov(y^​i​,yi​)/σ2 来验证．然而估计最优子集选择的 df(y^)\rm{df}(\hat{\mathbf y})df(y^​) 没有闭形式 (closed form)．

对于 LAR 和 lasso，会发生很奇怪的事情．这些技巧的自适应方式比最优集选择更加光滑，因此估计自由度会更加地难以驾驭．特别地，可以证明经过 kkk 步 LAR 过程，拟合向量的有效自由度恰巧是 kkk．对于 lasso，（改进的）LAR 过程经常需要多余 kkk 的步骤，因为可以删去预测变量．因此，定义有点不一样；对于 lasso，在任一小步 df(y^)\rm{df}(\hat{\mathbf y})df(y^​) 近似等于模型中预测变量的个数．然而这种近似在 lasso 路径中的任何地方都表现得很好，但是对于每个 kkk，它在包含 kkk 个预测变量的序列中最后一个模型表现得最好．关于 lasso 自由度详细的研究或许可以在 Zou et al. (2007)^{6 中找到．}

个人解读

关于有效自由度的定义式书中虽已经给出，但是如何理解呢，有的学者通过用它计算OLS和岭回归的自由度创造了感性理解的机会。

首先，下面等式中左右两部分是等价的变换
cov⁡(y^i,yi)=cov⁡(eiTy^,eiTy)=eiTcov⁡(y^,y)ei \operatorname{cov}\left(\hat{y}_{i}, y_{i}\right)=\operatorname{cov}\left(e_{i}^{T} \hat{y}, e_{i}^{T} y\right)=e_{i}^{T} \operatorname{cov}(\hat{y}, y) e_{i} cov(y^​i​,yi​)=cov(eiT​y^​,eiT​y)=eiT​cov(y^​,y)ei​
而对于最小二乘来说，预测变量的估计值y^=Xβ^ls=X(XTX)−1XTy\hat{y}=X \hat{\beta}^{\mathrm{ls}}=X\left(X^{T} X\right)^{-1} X^{T} yy^​=Xβ^​ls=X(XTX)−1XTy，所以上式可转变为
cov⁡(y^,y)=X(XTX)−1XTcov⁡(y,y)=σ2X(XTX)−1XT \operatorname{cov}(\hat{y}, y)=X\left(X^{T} X\right)^{-1} X^{T} \operatorname{cov}(y, y)=\sigma^{2} X\left(X^{T} X\right)^{-1} X^{T} cov(y^​,y)=X(XTX)−1XTcov(y,y)=σ2X(XTX)−1XT
上式最后两步是怎么得出来的，这是有关协方差的相关概念，cov⁡(y,y)=σ2I\operatorname{cov}(y, y)=\sigma^{2} Icov(y,y)=σ2I，因此
cov⁡(y^i,yi)=σ2eiTX(XTX)−1XTei=σ2(XTei)(XTX)−1(XTei) \operatorname{cov}\left(\hat{y}_{i}, y_{i}\right)=\sigma^{2} e_{i}^{T} X\left(X^{T} X\right)^{-1} X^{T} e_{i}=\sigma^{2}\left(X^{T} e_{i}\right)\left(X^{T} X\right)^{-1}\left(X^{T} e_{i}\right) cov(y^​i​,yi​)=σ2eiT​X(XTX)−1XTei​=σ2(XTei​)(XTX)−1(XTei​)
这里的XTei=xiX^Te_i=x_iXTei​=xi​，所以原式等于
df⁡(y^)=∑i=1NxiT(XTX)−1xi=∑i=1Ntrace⁡(xiT(XTX)−1xi)=∑i=1Ntrace⁡(xixiT(XTX)−1)=trace⁡((∑i=1NxixiT)(XTX)−1) \begin{aligned} \operatorname{df}(\hat{y}) &=\sum_{i=1}^{N} x_{i}^{T}\left(X^{T} X\right)^{-1} x_{i} \\ &=\sum_{i=1}^{N} \operatorname{trace}\left(x_{i}^{T}\left(X^{T} X\right)^{-1} x_{i}\right) \\ &=\sum_{i=1}^{N} \operatorname{trace}\left(x_{i} x_{i}^{T}\left(X^{T} X\right)^{-1}\right) \\ &=\operatorname{trace}\left(\left(\sum_{i=1}^{N} x_{i} x_{i}^{T}\right)\left(X^{T} X\right)^{-1}\right) \end{aligned} df(y^​)​=i=1∑N​xiT​(XTX)−1xi​=i=1∑N​trace(xiT​(XTX)−1xi​)=i=1∑N​trace(xi​xiT​(XTX)−1)=trace((i=1∑N​xi​xiT​)(XTX)−1)​
上式的第一步到第二步是因为最终结果是个标量，对标量取trace是其本身，第二步到第三步是因为trace(AB)=trace(BA)trace(AB)=trace(BA)trace(AB)=trace(BA)，最终结果里面的第一项又可以被写为
∑i=1NxixiT=[x1x2⋯xN][x1Tx2T⋮xNT]=XTX \sum_{i=1}^{N} x_{i} x_{i}^{T}=\left[\begin{array}{llll} {x_{1}} & {x_{2}} & {\cdots} & {x_{N}} \end{array}\right]\left[\begin{array}{c} {x_{1}^{T}} \\ {x_{2}^{T}} \\ {\vdots} \\ {x_{N}^{T}} \end{array}\right]=X^{T} X i=1∑N​xi​xiT​=[x1​​x2​​⋯​xN​​]⎣⎢⎢⎢⎡​x1T​x2T​⋮xNT​​⎦⎥⎥⎥⎤​=XTX
所以最终结果为
df⁡(y^)=trace⁡((XTX)(XTX)−1)=trace⁡(Ik×k)=k \operatorname{df}(\hat{y})=\operatorname{trace}\left(\left(X^{T} X\right)\left(X^{T} X\right)^{-1}\right)=\operatorname{trace}\left(I_{k \times k}\right)=k df(y^​)=trace((XTX)(XTX)−1)=trace(Ik×k​)=k
这就是最小二乘的有效自由度，那么岭回归的有效自由度推导如下，首先
y^=Xβ^ridge =X(XTX+λI)−1XTy \hat{y}=X \hat{\beta}^{\text {ridge }}=X\left(X^{T} X+\lambda I\right)^{-1} X^{T} \mathbf{y} y^​=Xβ^​ridge =X(XTX+λI)−1XTy
那么协方差就等于
cov⁡(y^,y)=X(XTX+λI)−1XTcov⁡(y,y)=σ2X(XTX+λI)−1XT \operatorname{cov}(\hat{y}, y)=X\left(X^{T} X+\lambda I\right)^{-1} X^{T} \operatorname{cov}(y, y)=\sigma^{2} X\left(X^{T} X+\lambda I\right)^{-1} X^{T} cov(y^​,y)=X(XTX+λI)−1XTcov(y,y)=σ2X(XTX+λI)−1XT
变成yiy_iyi​的形式就是
cov⁡(y^i,yi)=σ2(XTei)T(XTX+λI)−1(XTei)=σ2xiT(XTX+λI)−1xi \operatorname{cov}\left(\hat{y}_{i}, y_{i}\right)=\sigma^{2}\left(X^{T} e_{i}\right)^{T}\left(X^{T} X+\lambda I\right)^{-1}\left(X^{T} e_{i}\right)=\sigma^{2} x_{i}^{T}\left(X^{T} X+\lambda I\right)^{-1} x_{i} cov(y^​i​,yi​)=σ2(XTei​)T(XTX+λI)−1(XTei​)=σ2xiT​(XTX+λI)−1xi​
代入式(3.60)(3.60)(3.60)可得
df⁡(y^)=∑i=1Ntrace⁡(xixiT(XTX+λI)−1)=trace⁡(XTX(XTX+λI)−1)=trace⁡(X(XTX+λI)−1XT) \begin{aligned} \operatorname{df}(\hat{y}) &=\sum_{i=1}^{N} \operatorname{trace}\left(x_{i} x_{i}^{T}\left(X^{T} X+\lambda I\right)^{-1}\right) \\ &=\operatorname{trace}\left(X^{T} X\left(X^{T} X+\lambda I\right)^{-1}\right) \\ &=\operatorname{trace}\left(X\left(X^{T} X+\lambda I\right)^{-1} X^{T}\right) \end{aligned} df(y^​)​=i=1∑N​trace(xi​xiT​(XTX+λI)−1)=trace(XTX(XTX+λI)−1)=trace(X(XTX+λI)−1XT)​

Hoerl, A. E. and Kennard, R. (1970). Ridge regression: biased estimation for nonorthogonal problems, Technometrics 12: 55–67. ↩︎

Chen, S. S., Donoho, D. and Saunders, M. (1998). Atomic decomposition by basis pursuit, SIAM Journal on Scientific Computing 20(1): 33–61. ↩︎

Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society Series B. 67(2): 301–320. ↩︎

Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004). Least angle regression (with discussion), Annals of Statistics 32(2): 407–499. ↩︎

Osborne, M., Presnell, B. and Turlach, B. (2000a). A new approach to variable selection in least squares problems, IMA Journal of Numerical Analysis 20: 389–404. ↩︎

Zou, H., Hastie, T. and Tibshirani, R. (2007). On the degrees of freedom of the lasso, Annals of Statistics 35(5): 2173–2192. ↩︎

作者：Nstar-LDS

岭回归 lasso 学习笔记 拉格朗日乘子 学习 奇异值分解 回归