第三章 线性回归模型

3.1 导言

线性模型在 训练样本 (training cases) 很少、低信噪比或者稀疏数据的情况下可以很好地做预测，可以充分描述输入变量如何影响输出变量.

线性方法可以应用到输入变量的变换上.

3.2 线性回归模型和最小二乘法

线性回归模型有如下形式

\[ f(X)=\beta_0+\sum\limits_{j=1}^pX_j\beta_j \tag{3.1} \]

变量 X_j 来源有：

  • 定量的输入变量
  • 定量输入变量的变换, 比如对数变换, 平方根变换或者平方变换
  • 基函数展开, 比如 X2=X12,X_3=X_13, 产生一个多项式的表示
  • 定性输入变量水平 (level) 的数值或“虚拟”编码. 举个例子, 如果 G 是有 5 个水平的输入变量, 我们可能构造 Xj,j=1,,5 使得 Xj=I(G=j). 通过一系列独立于水平的常数, 整个 Xj 可以用来表示 G 的效应, 因为在 _{j=1}^5X_j_j 中, 其中一个 Xj 的系数为 1, 其它的都是 0.
  • 变量的交互影响, 举个例子, X3=X1X2.

无论 X_j 是哪个来源, 模型关于参数都是线性的.

最小二乘 (least squares), 我们取参数 =(_0,_1,,_p)^T 使得下式的残差平方和最小

\[\begin{align*} \text{RSS}(\beta)&=\sum\limits_{i=1}^N(y_i-f(x_i))^2\ &=\sum\limits_{i=1}^N(y_i-\beta_0-\sum\limits_{j=1}^px_{ij}\beta_j)^2\tag{3.2} \end{align*}\]s_{i=1}N(y_i-0-{j=1}px_{ij}_j)^2 \end{align*}

从统计学的观点来看, 如果训练观测值 (xi,yi) 是从总体独立随机抽取的, 则该标准是很合理的. 即使 xi’s 不是随机选取的, 如果在给定输入 xi 的条件下 yi 条件独立, 这个准则也是有效的.

·图 3.1 展示了在充满数据对 (X,Y) 的 ^{p+1} 维空间的最小二乘拟合的几何意义.

·式( 3.2 ) 对模型 (3.1) 的有效性没有作假设, 而是简单地找到对数据最好的线性拟合.

·无论数据是怎样产生的, 最小二乘拟合直观上看是满意的, 这个准则衡量了平均的拟合误差.

图 3.1 \(X\in R^2\) 中的最小二乘拟合. 我们寻找关于 \(X\) 并且使 \(Y\) 的残差最小的线性函数.

最小化式( 3.2 )

记 为 N(p+1) 的矩阵, 矩阵每一行为一个输入向量(在第一个位置为 1), 类似地令 为训练集里的 N 维输出向量.

残差平方和写成如下形式

\[ \text{RSS}(\beta)=(\mathbf{y}-\mathbf{X}\beta)^T(\mathbf{y}-\mathbf{X}\beta)\tag{3.3} \]

这是含 p+1 个参数的二次函数. 关于 求导有

\[ \begin{array}{ll} \dfrac{\partial \text{RSS}}{\partial \beta}&=-2\mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta)\\ \dfrac{\partial^2 \text{RSS}}{\partial \beta\partial \beta^T}&=2\mathbf{X}^T\mathbf{X} \end{array} \tag{3.4} \]

假设 列满秩, 因此 ^T 是正定的, 令一阶微分为 0, 即

\[ \mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta)=0\tag{3.5} \]

得唯一解

\[ \hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.6} \]

在输入向量 x_0 下的预测值由 (x_0)=(1:x_0)^T 给出；在训练输入值下的拟合值为

\[ \hat{y}=\mathbf{X}\hat{\beta}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.7} \]

其中, y_i=f(x_i). 在式( 3.7 ) 中出现的矩阵 “帽子”矩阵

图 3.2 含两个预测的最小二乘回归的 \(N\) 维几何图形. 输出向量 \(\mathbf{y}\) 正交投影在输入向量 \(\mathbf{x}_1\) 和\(\mathbf{x}_2\) 张成的超平面中. 投影向量 \(\hat{\mathbf{y}}\) 表示最小二乘法的预测值.

图 3.2是在 \(\mathbb{R}^N\) 中最小二乘估计的另一种几何表示. 记 \(\mathbf{X}\) 的列向量为 \(x_0,x_1,\ldots,x_p\), 其中 \(x_0 \equiv 1\). 这些向量张成了 \(\mathbb{R}^N\) 的子空间, 也被称作 \(\mathbf{X}\) 的列空间. 我们通过选择 \(\hat{\beta}\) 来最小化 \(\text{RSS}(\beta)=\lVert(\mathbf{y}-\mathbf{X}\beta)\rVert^2\), 则残差向量 \(\mathbf{y-\hat{y}}\) 正交于子空间. ( 3.5 ) 式描述了这种正交, 然后估计值 \(\hat{y}\) 因此可以看成是 \(\mathbf{y}\) 在子空间中的正交投影. 帽子矩阵 \(\mathbf{H}\) 计算了正交投影, 也就是指一类正交投影矩阵.

在出现 \(\mathbf{X}\) 非满秩的情形时，需要通过重编码或去除 \(\mathbf{X}\) 中冗余的问题.

假设观测值 \(y_i\) 不相关, 且有固定的方差 \(\sigma^2\), 并且 \(x_i\) 是固定的. 最小二乘法的参数估计的 协方差矩阵 为：

\[ \text{Var}(\hat{\beta})=(\mathbf{X}^T\mathbf{X})^{-1}\sigma^2\tag{3.8} \]

方差 \(\sigma^2\)通过下式来估计:

\[ \hat{\sigma}^2=\frac{1}{N-p-1}\sum\limits_{i=1}^N(y_i-\hat{y}_i)^2 \]

分母不是 \(N\), 是因为 \(N-p-1\) 使得 \(\hat{\sigma}^2\) 为无偏估计：\(E(\hat{\sigma}^2)=\sigma^2\)

假设式( 3.1 ) 是关于均值的正确模型；则 \(Y\) 的条件期望关于 \(X_1,X_2,\ldots,X_p\) 是线性的. 假设 \(Y\) 与其期望的偏差是可加的且是正态的. 因此 \[\begin{align*} Y&=\text{E}(Y\mid X_1,\ldots,X_p)+\varepsilon\\ &=\beta_0+\sum\limits_{j=1}^pX_j\beta_j+\varepsilon\tag{3.9} \end{align*}\] 其中误差 \(\epsilon\) 是期望值为 0 方差为 \(\sigma^2\) 的高斯随机变量, 记作 \(\varepsilon\sim N(0,\sigma^2)\)

由式( 3.9 ), 可以证明 \[ \hat{\beta}\sim N(\beta,(\mathbf{X}^T\mathbf{X})^{-1}\sigma^2) \tag{3.10} \] 这是一个有上述均值向量和协方差矩阵的多元正态分布. 同时有 \[ (N-p-1)\hat \sigma^2\sim \sigma^2\chi^2_{N-p-1}\tag{3.11} \] 是一个自由度为 \(N-p-1\) 的卡方分布. 另外, \(\hat{\beta}\) 和 \(\hat\sigma^2\) 是统计独立的. 我们利用这些分布性质得到假设检验以及对于参数 \(\beta_j\) 的置信区间

为检验系数 \(\beta_j\) 的假设, 构造标准化因数或者 Z-分数. \[ z_j=\dfrac{\hat{\beta}_j}{\hat{\sigma}\sqrt{v_j}}\tag{3.12} \] 当 \(z_j\) 的绝对值较大时会拒绝零假设. 需要同时检验多个系数的显著性时利用 \(F\) 统计量 \[ F=\dfrac{(\text{RSS}_0-\text{RSS}_1)/(p_1-p_0)}{\text{RSS}_1/(N-p_1-1)}\tag{3.13} \] 对于线性模型(比如多重回归), 可以用来进行 power analysis, 其中 u 和 v 分别是分子和分母的自由度, f2 是 有效大小. 效力函数可以帮助确定在指定显著性条件下实验所需要的样本量并评估实验设计的统计效力.

这里举一个实际例子来说明如何在线性回归模型中进行 power analysis.

假设我们要研究老板的领导风格对公司员工满意度的影响, 同时员工满意度还与薪酬有关, 其中薪酬用 3 个变量来衡量, 而老板领导风格用 4 个变量来衡量. 根据以往经验, 薪酬能够解释员工满意度方差的 30%, 现在问题是确定多少的调查个体, 使得显著度为 .05, 而且 power 为 .90. 求解代码为pwr.f2.test(u=3, f2=(.35-.30)/(1-.35), sig.level=0.05, power=0.90), 输出v=184.2426.

其中u和v分别被称为分子、分母自由度. 结合 式( 3.13 ) 式, 我们可以看出 \(p_1-p_0=u, N-p_1-1=v\), 并且这里 \(p_0=4, p_1=7\). 换个角度看, 其实式( 3.13 ) 就是用于 power analysis 的检验统计量, 原假设为可以将领导力水平的四个变量系数设为 0, 若 \(F\) 足够大, 则会拒绝原假设.

分离式( 3.10 ) 中的 \(\beta_j\) 得到 \(\beta_j\) 的置信水平为 \(1-2\alpha\) 的置信区间 \[ (\hat{\beta}_j-z^{1-\alpha}v_j^{1/2}\hat{\sigma},\hat{\beta}_j+z^{1-\alpha}v_j^{1/2}\hat{\sigma})\tag{3.14} \] 这里, \(z^{1-\alpha}\) 是正态分布的 \(1-\alpha\) 分位数. \[\begin{align*} z^{1-0.025}&=1.96,\\ z^{1-0.05}&=1.645,etc. \end{align*}\] 因此报告 \(\hat{\beta}\pm 2\cdot se(\hat{\beta})\) 的标准做法是约 \(95\%\) 的置信区间. 即便没有高斯误差的假设, 区间也近似正确, 因为当样本规模 \(N\sim \infty\) 时, 覆盖范围近似为 \(1-2\alpha\).

全体系数向量 \(\beta\) 的近似置信集, 记作 \[ F_{\beta} = \{\beta \mid ((\hat{\beta}-\beta)^T\mathbf{X^TX}(\hat{\beta}-\beta)\le \hat{\sigma}^2\chi^{2^{(1-\alpha)}}_{p+1} \}\tag{3.15} \] 其中, \(\chi_\ell^{2^{(1-\alpha)}}\) 是 \(\ell\) 个自由度的卡方分布的 \(1-\alpha\) 的分位数：举个例子, \(\chi^{2^{(1-0.05)}}_5=11.1,\chi^{2^{(1-0.1)}}_5=9.2\). 这个关于 \(\beta\) 的置信集导出关于真函数 \(f(x)=x^T\beta\) 对应的置信集, 记作 \(\{x^T\beta\mid \beta \in C_\beta\}\)

3.2.1 例子：前列腺癌

表 3.1 前列腺癌数据中预测变量的相关系数；表 3.2 前列腺癌数据的线性模型拟合. \(Z\) 分数为系数除以标准差 式( 3.12 ). 大体上, \(Z\) 分数绝对值大于 2 表示在 \(p=0.05\) 的水平下显著不为零.

测试数据的平均预测误差为 0.521. 相反地, 使用 lpsa 训练数据的平均值预测的测试误差为 1.057, 称作“基本误差阶”. 线性模型将降低基本误差阶.

3.2.2 Guass-Markov 定理

参数 \(\beta\) 的最小二乘估计在所有的线性无偏估计中有最小的方差. 考虑任意参数 \(\theta = a^T\beta\) 的线性组合, 举例: 预测值 \(f(x_0)=x_0^T\beta\) 是这种形式. \(a^T\beta\) 的最小二乘估计为 \[ \hat{\theta}=a^T\hat{\beta}=a^T\mathbf{(X^TX)^{-1}X^Ty}\tag{3.17} \] 考虑固定 \(\mathbf{X}\), 则上式是关于响应变量 \(\mathbf{y}\) 的线性函数 \(\mathbf{c}_0^T\mathbf{y}\). 如果我们假设线性模型是正确的, 则 \(a^T\hat{\beta}\) 是无偏的, 因为 \[\begin{align*} \text{E}(a^T\hat{\beta})&=E(a^T\mathbf{(X^TX)^{-1}X^Ty})\\ &=a^T\mathbf{(X^TX)^{-1}X^TX}\beta\\ &=a^T\beta\tag{3.18} \end{align*}\] 依据 Gauss-Markov 定理说如果 \(a^T\beta\) 存在其它无偏估计 \(\tilde{\theta}=\mathbf{c}^T\mathbf{y}\), 即 \(\text{E}(\mathbf{c}^T\mathbf{y})=a^T\beta\), 则 \[ Var(a^T\hat{\beta})\le Var(\mathbf{c}^T\mathbf{y})\tag{3.19} \] 用整个参数向量 \(\beta\) 来说明Gauss-Markov 定理.

考虑 \(\theta\) 的估计值 \(\tilde{\theta}\) 的均方误差 \[\begin{align*} \text{MSE}(\tilde{\theta})&=\text{E}(\tilde{\theta}-\theta)^2\\ &=\text{Var}(\tilde{\theta})+[\text{E}(\tilde{\theta})-\theta]^2\tag{3.20} \end{align*}\] 第一项为方差, 第二项为平方偏差.

Gauss-Markov 定理表明最小二乘估计在所有无偏线性估计中有最小的均方误差.

任何收缩或者将最小二乘的一些参数设为 0 的方法都可能导致有偏估计. 挑选一个合适的模型意味着要在偏差和方差之间创造某种良好的平衡. 均方误差与预测的正确性紧密相关. 考虑在输入 \(x_0\) 处的新的响应变量 \[ Y_0=f(x_0)+\varepsilon_0 \tag{3.21} \] 其估计量 \(\tilde{f}(x_0)=x_0^T\tilde{\beta}\) 的期望预测误差为 \[\begin{align*} \text{E}(Y_0-\tilde{f}(x_0))^2 &=\sigma^2+E(x_0^T\tilde{\beta}-f(x_0))^2\\ &= \sigma^2+MSE(\tilde{f}(x_0))\tag{3.22} \end{align*}\] 预测误差的估计值和均方误差只有常数值 \(\sigma^2\) 的差别, \(\sigma^2\) 表示了新观测 \(y_0\) 的方差.

3.2.3 从简单单变量回归到多重回归

多重线性回归模型指有 p > 1 个输入的线性模型 (3.1).

假设有一个没有截距的单变量模型 \[ Y=X\beta + \varepsilon \tag{3.23} \] 最小二乘估计和残差为 \[ \begin{aligned} \hat{\beta}&=\dfrac{\sum_1^Nx_iy_i}{\sum_1^Nx_i^2}\\ r_i &= y_i -x_i\hat{\beta} \end{aligned} \tag{3.24} \] 为了简便用向量表示, 我们令 \(\mathbf{y}=(y_1,\ldots,y_N)^T\), \(\mathbf{x}=(x_1,\ldots,x_N)^T\), 并且定义 \[\begin{align*} \langle\mathbf{x},\mathbf{y}\rangle &= \sum\limits_{i=1}^Nx_iy_i\\ &=\mathbf{x^Ty}\tag{3.25} \end{align*}\] \(\mathbf{x}\) 和 \(\mathbf{y}\) 之间的内积, 可以写成 \[ \begin{aligned} \hat{\beta}&=\dfrac{\langle \mathbf{x,y}\rangle}{\langle\mathbf{x,x} \rangle}\\ \mathbf{r}&=\mathbf{y}-\mathbf{x}\hat{\beta} \end{aligned} \tag{3.26} \]

内积表示是线性回归模型一般化到不同度量空间(包括概率空间)建议的方式.

这个简单的单变量回归提供了多重线性回归的框架.

假设输入变量 \(\mathbf{x}_1,\mathbf{x_2,\ldots,x_p}\)(数据矩阵 \(\mathbf{X}\) 的列)是正交的得到多重最小二乘估计 \(\hat{\beta}_j\) 等于 \(\langle \mathbf{x}_j,\mathbf{y}\rangle/\langle\mathbf{x}_j,\mathbf{x}_j\rangle\) ——单变量估计.

当输入变量为正交的, 它们对模型中其它的参数估计没有影响.

进行正交化, 假设我们有一个截距和单输入 \(\bf{x}\). 则 \(\bf{x}\) 的最小二乘系数有如下形式 \[ \hat{\beta}_1=\dfrac{\langle \mathbf{x}-\bar{x}\mathbf{1},\mathbf{y}\rangle}{\langle \mathbf{x}-\bar{x}\mathbf{1},\mathbf{x}-\bar{x}\mathbf{1}\rangle}\tag{3.27} \] 其中, \(\bar{x}=\sum_ix_i/N\),且 \(N\) 个元素全为 1 的向量 \(\mathbf{1}=\mathbf{x_0}\). 将式( 3.27 ) 的估计看成简单回归 式( 3.26 ) 的两次应用. 这两步是：

  1. 在 \(\bf{1}\) 上回归 \(\bf{x}\) 产生残差 \(\mathbf{z}=\mathbf{x}-\bar{x}\mathbf{1}\);
  2. 在残差 \(\bf{z}\) 上回归 \(\bf{y}\) 得到系数 \(\hat{\beta}_1\) 在这个过程中, “在 \(\bf{a}\) 上回归 \(\bf{b}\)”意思是 \(\bf{b}\) 在 \(\bf{a}\) 上的无截距的简单单变量回归, 产生系数 \(\hat{\gamma}=\langle\mathbf{a,b}\rangle/\langle\mathbf{a,a}\rangle\) 以及残差向量 \(\mathbf{b}-\hat{\gamma}\mathbf{a}\). 我们称 \(\bf{b}\) 由 \(\bf{a}\) 校正(adjusted), 或者关于 \(\bf{a}\) 正交化.

第一步对 \(\mathbf{x}\) 作关于 \(\mathbf{x}_0=\mathbf{1}\) 的正交化.

第二步是一个利用正交预测变量 \(\mathbf{1}\) 和 \(\mathbf{z}\) 简单的单变量回归. 图 3.4 展示了两个一般输入 \(\mathbf{x}_1\) 和 \(\mathbf{x}_2\) 的过程. 正交化不会改变由 \(\mathbf{x}_1\) 和 \(\mathbf{x}_2\) 张成的子空间, 它简单地产生一个正交基来表示子空间.

正交输入的最小二乘回归. 向量 \(\mathbf{x}_2\) 在向量 \(\mathbf{x}_1\) 上回归, 得到残差向量 \(\mathbf{z}\). \(\mathbf{y}\) 在\(\mathbf{z}\) 上的回归给出 \(\mathbf{x}_2\) 的系数. 把 \(\mathbf{y}\) 在 \(\mathbf{x}_1\) 和 \(\mathbf{z}\) 上的投影加起来给出了最小二乘拟合 \(\mathbf{\hat{y}}\).

推广到 \(p\) 个输入的情形, 如算法 3.1 所示. 计算得到的简单回归的系数实际上是多重回归的系数.

结果为 \[ \hat{\beta}_p=\dfrac{\langle \mathbf{z}_p,\mathbf{y} \rangle}{\langle \mathbf{z}_p,\mathbf{z}_p \rangle}\tag{3.28} \] 对第二步中的残差进行重排列, 我们可以看到每个 \(\mathbf x_j\) 是 \(\mathbf z_k,k\le j\) 的线性组合. 因为 \(\mathbf z_j\) 是正交的, 它们形成了 \(\mathbf{X}\) 列空间的基, 从而得到最小二乘在该子空间上投影为 \(\hat{\mathbf{y}}\). 因为 \(\mathbf z_p\)仅与 \(\mathbf x_p\) 有关(系数为 1), 则式( 3.28 ) 确实是 \(\mathbf{y}\) 在 \(\mathbf x_p\) 上多重回归的系数.

通过对 \(\mathbf x_j\) 重排列, 其中的任何一个都有可能成为最后一个, 然后得到一个类似的结果. 多重回归的第 \(j\) 个系数是 \(\mathbf{y}\) 在 \(\mathbf x_{j\cdot 012...(j-1)(j+1)...,p}\) 的单变量回归, \(\mathbf x_{j\cdot 012...(j-1)(j+1)...,p}\) 是 \(\mathbf x_j\) 在\(\mathbf x_0,\mathbf x_1,...,\mathbf x_{j-1},\mathbf{j+1},...,\mathbf{x}_p\) 回归后的残差向量：

多重回归系数 \(\hat\beta_j\) 表示 \(\mathbf x_j\) 经过 \(\mathbf x_0,\mathbf x_1,\ldots,\mathbf x_{j-1},\mathbf x_{j+1},\ldots,\mathbf x_p\) 调整后 \(\mathbf x_j\) 对 \(\mathbf y\) 额外的贡献

如果 \(\mathbf{x}_p\) 与某些 \(\mathbf{x}_k\) 高度相关, 残差向量 \(\mathbf{z}_p\) 会近似等于 0, 而且由 式( 3.28 ) 得到的系数 \(\hat{\beta}_p\) 会非常不稳定. 对于相关变量集中所有的变量都是正确的. 在这些情形下我们可能得到所有的 Z 分数(如表 3.2 所示)都很小. 由 式( 3.28 ) 我们也得到估计方差 式( 3.8 ) 的另一种方法 \[ \text{Var}(\hat{\beta}_p)=\dfrac{\sigma^2}{\langle \mathbf{z}_p,\mathbf{z}_p \rangle}=\dfrac{\sigma^2}{\Vert\mathbf{z}_p\Vert ^2}\tag{3.29} \] 换句话说, 我们估计 \(\hat{\beta}_p\) 的精度取决于残差向量 \(\mathbf{z}_p\) 的长度；它表示 \(\mathbf{x}_p\) 不能被其他 \(\mathbf{x}_k\) 解释的程度.

算法 3.1 被称作多重回归的 Gram-Schmidt 正交化, 也是一个计算估计的有用的数值策略. 从中不仅可以得到 \(\hat{\beta}_p\), 而且可以得到整个多重最小二乘拟合,

可以用矩阵形式来表示算法 3.1 的第二步： \[ \mathbf{X=Z\Gamma}\tag{3.30} \] 其中 \(\mathbf{Z}\) 是 \(\mathbf{z_j}\)(按顺序)作为列向量的矩阵, \(\mathbf{\Gamma}\) 是值为 \(\hat\gamma_{kj}\) 的上三角矩阵. 引入第 \(j\) 个对角元 \(D_{jj}=\Vert \mathbf{z}_j \Vert\) 的对角矩阵 \(\mathbf{D}\), 我们有 \[\begin{align*} \mathbf{X}&=\mathbf{ZD^{-1}D\Gamma}\\ &=\mathbf{QR}\tag{3.31} \end{align*}\] 称为矩阵 \(\mathbf{X}\) 的 QR 分解. 这里 \(\mathbf{Q}\) 为 \(N\times (p+1)\) 正交矩阵, \(\mathbf{Q^TQ=I}\), 并且 \(\mathbf{R}\) 为 \((p+1)\times (p+1)\) 上三角矩阵.

\(\mathbf{QR}\) 分解表示了 \(\mathbf{X}\) 列空间的一组方便的正交基. 举个例子, 很容易可以看到, 最小二乘的解由下式给出 \[ \hat{\beta}=\mathbf{R^{-1}Q^Ty}\tag{3.32} \]

\[ \hat{\mathbf{y}}=\mathbf{QQ^Ty}\tag{3.33} \]

等式( 3.32 ) 很容易求解, 因为 \(\mathbf{R}\) 为上三角

3.2.4 多重输出

假设有多重输出 \(Y_1,Y_2,\ldots,Y_K\), 通过输入变量 \(X_0,X_1,X_2,\ldots,X_p\) 去预测.

假设对于每一个输出变量有线性模型 \[ Y_k=\beta_{0k}+\sum\limits_{j=1}^pX_j\beta_{jk}+\varepsilon_k \tag{3.34} \]

\[ = f_k(X)+\varepsilon_k \tag{3.35} \]

有 \(N\) 个训练情形时将模型写成矩阵形式 \[ \mathbf{Y=XB+E}\tag{3.36} \] 这里 \(\mathbf{Y}\) 为 \(N\times K\) 的响应矩阵, \(ik\) 处值为 \(y_{ik}\), \(\mathbf{X}\) 为 \(N\times (p+1)\) 的输入矩阵, \(\mathbf{B}\) 为 \((p+1)\times K\) 的系数矩阵 \(\mathbf{E}\) 为 \(N\times K\) 的误差矩阵. 单变量损失函数 式( 3.2 ) 的直接推广为 \[ \text{RSS}(\mathbf{B})=\sum\limits_{k=1}^K\sum\limits_{i=1}^N(y_{ik}-f_k(x_i))^2\tag{3.37} \]

\[ = \text{tr}[\mathbf{(Y-XB)^T(Y-XB)}]\tag{3.38} \]

最小二乘估计有和前面一样的形式 \[ \mathbf{\hat{B}=(X^TX)^{-1}X^TY}\tag{3.39} \] 因此第 \(k\) 个输出的系数恰恰是 \(\mathbf{y}_k\) 在 \(\mathbf{x}_0,\mathbf{x}_1,\ldots,\mathbf{x}_p\) 上回归的最小二乘估计. 多重输出不影响其他的最小二乘估计.

如果式( 3.34 ) 中的误差 \(\varepsilon = (\varepsilon_1,\ldots,\varepsilon_K)\) 相关, 则似乎更恰当的方式是修正式( 3.37 ) 来改成多重变量版本. 特别地, 假设 \(\text{Cov}(\varepsilon)=\Sigma\), 则 多重变量加权准则 为 \[ \text{RSS}(\mathbf{B;\Sigma})=\sum\limits_{i=1}^N(y_i-f(x_i))^T\Sigma^{-1}(y_i-f(x_i))\tag{3.40} \] 这可以由 多变量高斯定理 自然得出. 高斯分布可从多个角度应用于不同场景下.

3.3 选择预测变量的子集

经常不满足最小二乘估计的原因:

  • 一、预测的精确性
  • 二、可解释性

一些线性回归选择变量子集的方法属于 模型选择.

子集选择意味着我们只保留变量的一个子集, 并除去模型中的剩余部分. 最小二乘回归用来预测保留下的输入变量的系数. 这里有一系列不同的选择子集的策略.

3.3.1 最优集的选择

对于每个 \(k\in \{0,1,2,\ldots,p\}\), 最优子集回归要找出规模为 \(k\) 的子集中残差平方和 (3.2) 最小的子集.

例如前列腺癌例子中, 规模为 2 的最优子集不需要包含规模为 1 最优子集中的变量. 图中的红色下边界最优子集曲线必然下降, 所以不能用来选择子集的规模 \(k\). 怎样选择 \(k\) 涉及偏差和方差之间的平衡, 以及各种准则.

选择最小的模型使得预测误差期望值的估计最小.

这些方法是使用训练数据去得到区别于复杂度和由单参数编码的模型序列.

3.3.2 向前和向后逐步选择

向前逐步选择从截距开始, 然后向模型中依次添加最大程度提升拟合效果的预测变量.

类似最优子集回归, 向前逐步产生由 \(k\) 索引的模型序列, \(k\) 为子集规模, 也是必须要确定的值.

向前逐步选择是 贪心算法, 产生一个嵌套的模型序列.

向前逐步选择更优的原因：

  • 计算 ；对于大的 \(p\) 值, 我们不能计算最优子集序列, 但是可以计算向前逐步序列
  • 统计 ；在每个规模下选择最优子集需要在方差上付出代价；向前逐步是一种有更多约束的搜索, 而且将会有更低的方差.

向后逐步选择 从整个模型开始, 逐步删掉对拟合影响最低的预测变量. 要删掉的候选变量是 Z 分数最低的变量. 向后只能用于 \(N>p\) 时, 而向前逐步都可以使用.

图 3.6 展示了一个用于比较最优子集回归和简单的向前向后选择的小型模拟研究的结果. 它们非常相似

在前列腺癌例子中, 最优子集、向前和向后选择都给出了完全相同的变量序列.

3.3.3 向前逐渐 (Forward-Stagewise) 回归

向前逐渐回归比向前逐步回归有更多限制.

类似向前逐步回归的是, 由等于 \(\bar{y}\) 的截距开始, 中心化后的预测变量系数都初始化为 0. 每一步算法找出与当前残差最相关的变量. 然后计算所选择变量的残差的简单线性回归系数, 并且添加到该变量的当前系数. 这个过程一直继续直到没有变量与残差有相关性——比如, 当 \(N>p\) 时的最小二乘拟合.

不同于向前逐步回归的是, 当一个变量添加到模型中其他的变量不需要调整. 结果导致, 向前逐渐回归需要用多于 \(p\) 步达到最小二乘拟合, 也就是慢拟合，适用于高维问题.

3.3.4 前列腺癌例子(继续)

表 3.3 展示了一系列不同的选择和收缩方法的系数.

表3.3. 对前列腺数据应用不同子集和收缩方法估计的系数和测试误差结果. 空白项对应于省略的变量.

图 3.7 展示了估计的预测误差曲线.

3.4 收缩的方法

收缩方法因为更加连续所以不会受高易变性太大的影响.

3.4.1 岭回归

岭回归 根据回归系数的大小加上惩罚因子对它们进行收缩. 岭回归的系数使得带惩罚的残差平方和最小 \[ \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} \] 这里\(\lambda\ge 0\)是控制收缩程度的参数：\(\lambda\)值越大, 收缩的程度越大. 每个系数都向零收缩. 通过参数的平方和来惩罚的想法也用在了神经网络, 也就是权重衰减 .

岭回归问题可以等价地写成 \[ \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} \] 上式用参数显式表达了对回归参数大小的约束.

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

对输入按比例进行缩放时, 岭回归的解不相等, 因此求解式( 3.41 ) 前我们需要对输入进行标准化.

将 式( 3.41 ) 的准则写成矩阵形式 \[ \text{RSS}(\lambda)=(\mathbf{y}-\mathbf{X}\beta)^T(\mathbf{y}-\mathbf{X}\beta)+\lambda\beta^T\beta \tag{3.43} \] 可以简单地看出岭回归的解为 \[ \hat{\beta}^{ridge}=(\mathbf{X^TX}+\lambda \mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.44} \]

图3.8展示了前列腺癌例子的岭回归系数估计, 绘制成关于 \(\text{df}(\lambda)\) 的函数图, \(\text{df}(\lambda)\) 为由惩罚 \(\lambda\) 得到的有效自由度.

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

岭回归估计是后验分布的众数；因分布为高斯分布, 则也是后验分布的均值.

中心化输入矩阵 \(\mathbf{X}\) 的 奇异值分解 在许多统计方法分析中非常有用. \(N\times p\) 阶矩阵 \(\mathbf{X}\) 的 SVD 分解有如下形式 \[ \mathbf{X=UDV^T}\tag{3.45} \] 这里 \(\mathbf{U}\) 和 \(\mathbf{V}\) 分别是 \(N\times p\) 和 \(p\times p\) 的正交矩阵, \(\mathbf{U}\)的列张成 \(X\) 的列空间, \(\mathbf{V}\) 的列张成 \(X\) 的行空间. \(\mathbf{D}\) 为 \(p\times p\) 的对角矩阵, 对角元 \(d_1\ge d_2 \ge \cdots \ge d_p \ge 0\) 称作 \(\mathbf{X}\) 的奇异值. 如果一个或多个 \(d_j=0\), 则 \(\mathbf{X}\) 为奇异的.

\[ f(x,y)=x^TAy,\qquad A\in \mathbb{R}^{n\times m} \]

出发, 通过引入线性变换

\[ x=U\xi,\qquad y=V\eta \]

将双线性函数变为

\[ f(x,y)=\xi^TS\eta \]

其中

\[ S=U^TAV\, \]

若选择 \(U\) 和 \(V\) 为正交矩阵, 则他们的选择各存在 \(n^2-n\) 个自由度. 他提出利用这些自由度使矩阵 \(S\) 的非对角元为0, 即矩阵\(S=\Sigma=\text{diag}(\sigma_1,\sigma_2,\ldots,\sigma_n)\)为对角矩阵. 则

\[ A=U\Sigma V^T \]

任意复长方矩阵奇异值分解定理称为 Autonee-Eckart-Young 定理, 即

令 \(A\in \mathbb{R}^{m\times n}\)(或\(C^{m\times n}\)),则存在正交(或酉)矩阵 \(U\in \mathbb{R}^{m\times m}\)(或 \(C^{m\times m}\))和 \(V\in \mathbb{R}^{n\times n}\)(或\(C^{n\times n}\))使得

\[ A=U\Sigma V^T(or\quad U\Sigma V^H) \]

式中

\[ \Sigma= \left[ \begin{array}{cc} \Sigma_1&O\\ O&O \end{array} \right] \]

且 \(\Sigma_1=diag(\sigma_1,\sigma_2,\ldots,\sigma_r)\) , 其对角元素按照顺序

\[ \sigma_1> \sigma_2\cdots\ge\sigma_r>0,\qquad r=rank(A) \]

排列.

下图(来自维基百科)形象地展示了 SVD 的四种不同形式,

  • Full SVD:
  • Thin SVD: 只保留 \(U\) 中对应 \(V^T\) 中行向量的 \(n\) 个列向量
  • Compact SVD: 去掉零奇异值对应的行和列
  • Truncated SVD: 保留前 \(t\) 个最大奇异值对应的行和列

利用奇异值分解,把最小二乘拟合向量写成 \[\begin{align*} \mathbf{X}\hat{\beta}^{ls}&=\mathbf{X(X^TX)^{-1}X^Ty}\\ &=\mathbf{UU^Ty}\tag{3.46} \end{align*}\] \(\mathbf{U}^T\mathbf y\) 是 \(\mathbf{y}\) 正交基 \(\mathbf{U}\) 下的坐标.

note “weiya 注：Recall”

\[ \hat{\beta}=\mathbf{R^{-1}Q^Ty}\tag{3.32} \]

\[ \hat{\mathbf{y}}=\mathbf{QQ^Ty}\tag{3.33} \]

\(\mathbf{Q}\) 和 \(\mathbf{U}\) 是 \(\mathbf{X}\) 列空间的两个不同的正交基.

现在岭回归的解为 \[\begin{align*} \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{align*}\] 其中 \(\mathbf{u}_j\) 是 \(\mathbf{U}\) 的列向量. 注意到因为 \(\lambda\ge0\), 我们有 \(d_j^2/(d^2_j+\lambda)\le1\). 类似线性回归, 岭回归计算 \(\mathbf{y}\) 关于正规基 \(\mathbf{U}\) 的坐标. 通过因子 \(d^2_j/(d^2_j+\lambda)\) 来收缩这些坐标. 这意味着更小的 \(d_j^2\) 会在更大程度上收缩基向量的坐标. \[ \mathbf{X^T X = VD^2V^T} \tag{3.48} \] 上式是 \(\mathbf{X^TX}\)(当忽略因子 \(N\) 时, 也是 \(S\))的 特征值分解 (eigen decomposition). 特征向量 \(v_j\)(\(\mathbf{V}\) 的列向量)也称作 \(\mathbf{X}\) 的 主成分 (principal components)(或 Karhunen-Loeve)方向. 第一主成分方向 \(v_1\) 有下面性质：\(\mathbf{z}_1=\mathbf{X}v_1\) 在所有 \(\mathbf{X}\) 列的标准化线性组合中有最大的样本方差. 样本方差很容易看出来是 \[ \text{Var}(\mathbf{z}_1)=\text{Var}(\mathbf{X}v_1)=\dfrac{d_1^2}{N}\tag{3.49} \] 事实上 \(\mathbf{z}_1=\mathbf{X}v_1=\mathbf{u}_1d_1\). 导出变量 \(\mathbf{z_1}\) 称作 \(\mathbf{X}\) 的第一主成分, 因此 \(\mathbf{u_1}\) 是标准化的第一主成分. 后面的主成分 \(z_j\) 在与前一个保持正交的前提下有最大的方差 \(d_j^2/N\). 所以, 最后一个主成分有最小的方差. 因此越小的奇异值 \(d_j\) 对应 \(\mathbf{X}\) 列空间中方差越小的方向, 并且岭回归在这些方向上收缩得最厉害.

图 3.9 展示了两个维度下部分数据点的主成分. 如果在这个区域(\(Y\) 轴垂直纸面)内拟合线性曲面, 数据的结构形态使得确定梯度时长方向会比短方向更精确. 岭回归防止在短方向上估计梯度可能存在的高方差.

图 3.7 预测误差估计值关于 \(\text{df}(\lambda)\) 的曲线 \[\begin{align*} \text{df}(\lambda)&=\text{tr}[\mathbf{X}(\mathbf{X^TX}+\lambda\mathbf{I})^{-1}\mathbf{X}^T]\\ &=\text{tr}(\mathbf{H}_{\lambda})\\ &=\sum\limits_{j=1}^p\dfrac{d_j^2}{d_j^2+\lambda}\tag{3.50} \end{align*}\] 上面 \(\lambda\) 的单调递减函数是岭回归拟合的 有效自由度. 表 3.3 表明岭回归将全最小二乘估计的测试误差降低了一小部分.

3.4.2 Lasso

lasso 也是个收缩方法. lasso 估计定义如下 \[\begin{align*} \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{align*}\] 正如在岭回归中一样, 我们可以通过标准化预测变量来对常数 \(\beta_0\) 再参量化；\(\hat{\beta}_0\) 的解为 \(\bar{y}\), 并且后面我们拟合无截距的模型.

可以把 lasso 问题等价地写成 拉格朗日形式 \[ \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可以解决高维数据的稀疏性问题.可以把不重要变量的系数压缩为0.既实现了较为准确的参数估计也实现了降维. 图 3.10 显示了当惩罚参数 \(s=t/\sum_1^p\vert\hat{\beta}_j\vert\) 不同时的 lasso 系数. 当 \(s=1.0\) 时为最小二乘估计；当 \(s\rightarrow 0\) 时下降为 0. 该下降不总是严格单调的.

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

有约束的线性回归模型的三种方法：子集选择、岭回归和 lasso.

如表3.4 每种方法对最小二乘估计 \(\hat{\beta}_j\) 应用简单的变换.

岭回归做等比例的收缩.

lasso 通过常数因子 \(\lambda\) 变换每个系数, 在 0 处截去.

最优子集选择删掉所有系数小于第 \(M\) 个大系数的变量.

当只有两个参数时图 3.11 描绘了 lasso(左)和岭回归(右). 残差平方和为椭圆形的等高线, 以全最小二乘估计为中心. 岭回归的约束区域为圆盘 \(\beta_1^2+\beta_2^2\le t\), lasso 的约束区域为菱形\(\vert\beta_1\vert+\vert\beta_2\vert\le t\). 两种方式都寻找当椭圆等高线达到约束区域的第一个点. 与圆盘不同, 菱形 有角；如果解出现在角上, 则有一个参数 \(\beta_j\) 等于 0. 当 \(p > 2\), 菱形变成了 偏菱形.

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

把岭回归和 lasso 一般化, 并且可以看成是贝叶斯估计. 考虑下面准则 \[ \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} \] 其中 \(q\ge 0\). 图 3.12 显示了两个输入情形下常数值 \(\sum_j\vert\beta_j\vert^q\) 的等高线.

lasso、岭回归和最优子集选择是有着不同先验分布的贝叶斯估计. 在贝叶斯估计中使用后验分布的均值更加常见. 岭回归同样是后验分布的均值, 但是 lasso 和最优子集选择不是. \[ \lambda \sum\limits_{j=1}^p(\alpha\beta_j^2+(1-\alpha)\vert\beta_j\vert)\tag{3.54} \] 这是一种岭回归和 lasso之间的平衡. 图 3.13 比较了 \(q=1.2\) 下的 \(L_q\) 惩罚以及 \(\alpha=0.2\) 的弹性网惩罚；很难从肉眼来观察出差异. 弹性网像 lasso 一样选择变量, 同时像岭回归一样收缩相关变量的系数. 同时考虑了 \(L_q\) 惩罚的计算优势.

3.4.4 最小角回归

最小角回归 可以看成是一种向前逐步回归的一个版本.

最小角回归和向前逐步回归差不多的策略, 但是仅仅加入一个变量应有的程度.

算法 3.2 最小角回归

  1. 对预测变量进行标准化处理得到零均值和单位范数. 以残差向量 \(\mathbf{r=y-\bar{y}},\beta_1,\ldots,\beta_p=0\) 开始.
  2. 找出与 \(\mathbf{r}\) 最相关的预测变量 \(\mathbf x_j\)
  3. 从 0 开始移动 \(\beta_j\) 一直到最小二乘系数 \(\langle\mathbf x_j, \mathbf r\rangle\), 直到存在其它的预测变量 \(\mathbf x_k\) 使得其与当前残差的相关性等于 \(\mathbf x_j\) 与当前残差的相关性.
  4. 在由当前残差在 \((\mathbf x_j,\mathbf x_k)\) 上的联合最小二乘系数方向上移动 \(\beta_j\) 和 \(\beta_k\), 直到存在其它的预测变量 \(x_l\) 与当前残差的相关性和当前残差与 \((\mathbf x_j,\mathbf x_k)\) 的相关性相等.
  5. 按这种方式继续直到所有的 \(p\) 个预测变量加入到模型中. 经过 \(\min(N-1, p)\) 步, 我们达到了全最小二乘的解.

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

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

由构造知 LAR 的系数以一种分段线性的方式进行改变. 图 3.15(左图)显示了 LAR 系数曲线作为 \(L_1\) 弧长的函数曲线.

图 3.15 的右图展示了对同样数据的 lasso 系数曲线. 几乎与左图相同, 当绿色曲线通过 0 时首次出现不同.

例如前列腺癌数据, LAR 系数曲线显示与图 3.10 的 lasso 曲线相同, 该曲线从不经过 0. 这些观测值促使对 LAR 算法进行简单修改, 给出了整个 lasso 路径, 它同样也是分段线性的.

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

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

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

图 3.16 将 LAR 和 lasso 与向前逐步和向前逐渐回归进行比较 . LAR 和 lasso 的行为与向前逐渐回归相似. 增长的向前逐渐回归与 LAR 和 lasso 类似

LAR 和 Lasso 自由度公式

1、首先考虑采用 \(k\) 个特征的子集的线性回归.

2、然后在该拟合模型中的自由度定义为\(k\).

我们需要一个对于自适应拟合模型的有效自由度的一般定义. 我们定义拟合向量 \(\hat{\mathbf y}=(\hat y_1,\hat y_2,\ldots,\hat y_N)\) 的自由度为 \[ \text{df}(\hat{\mathbf y})=\frac{1}{\sigma^2}\sum\limits_{i=1}^N\text{Cov}(\hat y_i,y_i)\tag{3.60} \] 这里 \(\text{Cov}(\hat y_i,y_i)\) 指的是预测值 \(\hat y_i\) 和其对应的输出向量 \(y_i\) 之间的协方差. 直观上看有意义：当拟合数据越困难, 协方差会越大, 从而 \(\text{df}(\hat{\mathbf y})\) 越大.

现在对于有 \(k\) 个固定预测变量的线性回归模型, 可以简单地证明 \(\text{df}(\hat{\mathbf y})=k\). 同样地, 对于岭回归, 这一定义导出表达式(3.50)的 闭型解：\(\text{df}(\hat{\mathbf{y}})=\text{tr}(\mathbf S_\lambda)\).

3.5 运用派生输入方向的方法

3.5.1 主成分回归

主成分回归构造派生的输入列 \(\mathbf z_m=\mathbf Xv_m\), 然后在 \(\mathbf z_1,\mathbf z_2,\ldots,\mathbf z_M,\; M\le p\) 上回归 \(\mathbf y\). 因为 \(\mathbf z_m\) 是正交的, 则这个回归只是单变量回归的和 \[ \hat{\mathbf y}^{pcr}_{(M)}=\bar y\mathbf 1+\sum\limits_{m=1}^M\hat{\theta}_m\mathbf z_m\tag{3.61} \] 其中, \(\hat\theta_m=\langle \mathbf z_m,\mathbf y\rangle/\langle\mathbf z_m,\mathbf z_m\rangle\). 因为每个 \(\mathbf z_m\) 是原输入变量 \(\mathbf x_j\) 的线性组合, 我们可以将解 式( 3.61 ) 表达成关于 \(\mathbf x_j\) 的系数: \[ \hat\beta^{pcr}(M)=\sum\limits_{m=1}^M\hat\theta_mv_m\tag{3.62} \]

主成分回归与岭回归非常相似：都是通过输入矩阵的主成分来操作的. 岭回归对主成分系数进行了收缩, 收缩更多地依赖对应特征值的大小；主成分回归丢掉 \(p-M\) 个最小的特征值分量. 如图 3.17

3.5.2 偏最小二乘

这个方法用于回归的输入变量的线性组合, 但是与主成分回归不同的是它采用 \(\mathbf y\)(除了 \(\mathbf X\))来构造. 和主成分回归相同的是, 偏最小二乘 (PLS) 也不是尺度不变 的 > note “weiya 注：” > 在 \(\mathbf a\) 上回归 \(\mathbf b\)(或者称作 \(\mathbf b\) 在 \(\mathbf a\) 上回归)指的是 > \(\mathbf b\) 在 \(\mathbf a\) 上的无截距的简单单变量回归, 回归系数为

\[ \hat \gamma = \dfrac{\langle \mathbf a,\mathbf b\rangle}{\langle \mathbf a,\mathbf a\rangle} \]

第 \(m\) 个主成分方向 \(v_m\) 是下面问题的解： \[ \begin{aligned} \max_{\alpha} &\text{ Var}(\mathbf X\alpha)\\ \text{ subject to } &\Vert \alpha\Vert=1,\alpha^T\mathbf Sv_\ell=0,\;\ell=1,\ldots,m-1 \end{aligned} \tag{3.63} \] \(\mathbf S\) 为 \(\mathbf x_j\) 的样本协方差矩阵. \(\alpha^T\mathbf Sv_\ell=0\) 保证了 \(\mathbf z_m=\mathbf X\alpha\) 与之前所有的线性组合 \(\mathbf z_\ell=\mathbf v_\ell\) 都不相关. 第 \(m\) 个 PLS 方向 \(\hat \varphi_m\) 是下面的解： \[ \begin{aligned} \max_\alpha & \text{ Corr}^2(\mathbf y,\mathbf X\alpha)\text{Var}(\mathbf X\alpha)\\ \text{ subject to } & \Vert\alpha\Vert=1,\alpha^T\mathbf S\hat \varphi_\ell=0,\ell=1,\ldots,m-1 \end{aligned} \tag{3.64} \] 方差项趋向于占主导地位, 偏最小二乘表现得很像岭回归和主成分回归. 偏最小二乘适用于处理预测变量之间存在多重共线性或预测变量数量多于观测值的情况. # 3.6 讨论：选择和收缩方法的比较

理解上面描述的不同方法之间的关系. 考虑相关系数为 \(\rho\) 的两个相关输入变量\(\mathbf X_1,\mathbf X_2\). 我们假设实际的回归系数为 \(\beta_1=4,\beta_2=2\).

上图显示了不同方法下当它们惩罚参数改变时的系数曲线.

总结到对于最小化预测误差, 岭回归一般比变量子集选择、主成分回归和偏最小二乘更好.

PLS, PCR 以及岭回归趋向于表现一致. 岭回归可能会更好, 因为它收缩得很光滑, 不像离散步骤中一样. Lasso 介于岭回归和最优子集回归中间, 并且有两者的部分性质.

3.7 多输出的收缩和选择

多重输出线性模型的最小二乘估计是关于每个输出的最小二乘估计.

note “weiya 注：”

\[ \hat{\mathbf \eta}^{ridge}=(\mathbf{X^TX}+\lambda \mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}\tag{3.44} \]

为了在多重输出情形中应用选择和收缩的方法, 可以对每个输出变量单独地应用单变量的技巧, 或对全部的输出变量同时进行.

例如, 假设在这些输出变量中, 有 \[ Y_k = f(X) + \varepsilon_k\tag{3.65} \]

\[ Y_\ell = f(X) + \varepsilon_\ell\,;\tag{3.66} \]

也就是说, 式( 3.65 ) 和 式( 3.66 ) 在它们的模型中共享相同的结构 \(f(X)\). 这种情形下很明显应该合并 \(Y_k\) 和 \(Y_\ell\) 来估计共同的 \(f\).

合并响应变量是 典则相关分析 的核心, 这是一种为多元输出情形提出的数据降维的技巧. 类似 PCA, CCA 寻找 \(x_j\) 的一列不相关的线性组合 \(X v_m, m=1,\ldots,M\), 以及对应的响应变量 \(y_k\) 的不相关的线性组合 \(Y u_m\), 使得相关系数 \[ \text{Corr}^2(Y u_m,X v_m)\tag{3.67} \] 降秩回归 采用显式地合并信息的回归模型来 正式化 这个方法. 给定误差协方差 \(\text{Cov}(\varepsilon)=\mathbf{\Sigma}\) , 求解下列带约束的多元回归问题： \[ \hat{\mathbf{B}}^{rr}(m) = \operatorname*{argmin}_{\text{rank}(\mathbf{B})=m} \sum_{i=1}^N(y_i-\mathbf{B}^Tx_i)^T\mathbf{\Sigma}^{-1}(y_i-\mathbf{B}^Tx_i). \tag{3.68} \] 将 \(\mathbf{\Sigma}\) 用估计值 \(Y^TY/N\) 替换, 可以证明(练习 3.21)解由 \(Y\) 和 \(X\) 的 CCA 给出： \[ \hat{\mathbf{B}}^{rr}(m)=\hat{\mathbf{B}} \mathbf{U}_m \mathbf{U}_m^- \tag{3.69} \] 其中 \(\mathbf{U}_m\) 是 \(\mathbf{U}\) 中由前 \(m\) 列构成的 \(K \times m\) 的子矩阵, \(\mathbf{U}\) 是 \(K \times M\) 的左典则向量 \(u_1,u_2,\ldots,u_M\) 构成的矩阵. \(\mathbf{U}_m^-\) 是其广义逆.

将解改写成 \[ \hat{\mathbf{B}}^{rr}(M) = (X^TX)^{-1}X^T(Y\mathbf{U}_m)\mathbf{U}_m^- \tag{3.70} \] 降秩回归在合并的响应矩阵 \(Y{\mathbf{U}}_m\) 上进行线性回归, 然后将系数映射回原来的响应变量空间中(因此能拟合得很好). 降秩拟合由下式给出： \[\begin{align*} \hat{Y}^{rr}(m) &= X(X^TX)^{-1}X^TY\mathbf{U}_m \mathbf{U}_m^-\\ &=HYP_m \tag{3.71} \end{align*}\] 其中 \(\mathbf{H}\) 是一般的线性回归映射算子, 而 \(\mathbf{P}_m\) 是秩为 \(m\) 的 CCA 响应变量投影算子. 尽管 \(\mathbf{\Sigma}\) 更好的估计为 \((Y-X\hat{\mathbf{B}})^T(Y-X\hat{\mathbf{B}})/(N-pK)\), 但是可以证明解是一样的.

降秩回归通过截断 CCA 从响应变量中借来了优势变量. \[ \hat{\mathbf{B}}^{c+w} = \hat{\mathbf{B}}{\mathbf{U}}\mathbf \Lambda{\mathbf{U}}^{-1}\,,\tag{3.72} \] 其中 \(\mathbf \Lambda\) 是对角收缩矩阵. 基于总体设定中的最优预测, 他们证明了 \(\mathbf \Lambda\) 的对角元为 \[ \lambda_m = \frac{c_m^2}{c_m^2+\frac pN(1-c_m^2)}\,,\; m = 1,\ldots,M \tag{3.73} \] 其中 \(c_m\) 是第 \(m\) 个典则相关系数. 注意到随着输入变量个数与样本大小的比率 \(p/N\) 变小, 收缩因子趋向于 1. 基于训练数据和交互验证提出了修改版本的 \(\mathbf \Lambda\), 一般形式是相同的. 这里拟合响应形式为 \[ \hat{Y}^{c+w} = HYS^{c+w} \tag{3.74} \] 这里 \(\S^{c+w} = {\mathbf{U}}\mathbf \Lambda{\mathbf{U}}^{-1}\) 是响应变量收缩算子.

Breiman and Friedman (1997)[^3] 还建议同时在 \(Y\) 和 \(X\) 的空间中进行收缩. 这导出混合收缩模型： \[ \hat{Y}^{\text{ridge},c+w}=A_\lambda YS^{c+w} \tag{3.75} \] 其中 \(A_\lambda = X(X^TX+\lambda I)^{-1}X^T\) 是岭回归收缩算子.

3.8 Lasso 和相关路径算法的补充

\(L_1\) 正则促进了信号处理领域的 压缩传感 的发展.

3.8.1 增长的向前逐渐回归 (Incremental Forward Stagewise Regression)

它通过重复更新与当前残差最相关的变量的系数(乘以一个小量 \(\epsilon\))得到系数曲线. 算法 3.4 给出了具体的细节. 图 3.19(左图)展示了前列腺癌数据中步长 \(\epsilon=0.01\) 的过程. 如果 \(\delta_j=\langle \mathbf x_j,\mathbf r\rangle\)(残差在第 \(j\) 个预测变量的最小二乘系数), 这就是 向前逐渐过程 (FS).

研究小的 \(\epsilon\) 值 令 \(\epsilon\rightarrow 0\) 则得到图 3.19 的右图, 这个极限过程为 无穷小的向前逐渐回归 或者 \(FS_0\) .

LAR 在这些连结预测变量中的最小二乘拟合可以导致系数向相反的方向移动到它们的相关系数, 对 LAR 算法的修正实现了 \(FS_0\)：

这个修正相当于一个非负的最小二乘拟合, 保持系数的符号与相关系数的符号一致. 可以证明它实现了对于最大相关性的连结变量的无限小”更新”的最优平衡.

结果来看, 如果 LAR 图象是单调不减或者单调不增, 如 3.19 所示, 则 LAR, lasso, 以及 \(FS_0\) 这三种算法给出了相同的图象. 如果图象不是单调的但是不穿过 0, 则 LAR 和 lasso 是一样的.

\(FS_0\) 比 lasso 更加的复杂；\(FS_0\) 系数曲线是微分方程的一个解. 尽管 lasso 在降低系数向量 \(\beta\) 的 \(L_1\) 范数的单位残差平方和增长方面实现了最优化, 但 \(FS_0\) 在沿着系数路径的 \(L_1\) 弧长的单位增长是最优的. 它的系数曲线不会经常改变方向.

\(FS_0\) 比 lasso的约束更强, 见图 16.3 . \(FS_0\) 可能在 \(p>>N\) 情形下很有用, 它的系数曲线会更加的光滑, 比 lasso 有更小的方差.

3.8.2 分段线性路径算法

最小角回归过程探索了 lasso 解的路径分段线性的本质. 这导出了其他正则化问题类似的“路径算法”. 假设求解 \[ \hat \beta(\lambda)=\text{argmin}_\beta[R(\beta)+\lambda J(\beta)]\tag{3.76} \]

\[ F(\beta)=\sum\limits_{i=1}^NL(y_i,\beta_0+\sum_{j=1}^px_{ij}\beta_j)\tag{3.77} \]

其中损失函数 \(L\) 和惩罚函数 \(J\) 都是凸函数. 则下面是解的路径 \(\hat\beta(\lambda)\) 为分段线性的充分条件

  1. \(R\) 作为 \(\beta\) 的函数是二次的或者是分段二次
  2. \(J\) 关于 \(\beta\) 分段线性

note “weiya 注：Huber Loss & Hinge Loss” Huber loss 损失函数为：

\[ \rho(t;\lambda) = \begin{cases} \lambda \vert t\vert - \lambda^2/2 & \text{ if }\vert t\vert >\lambda\\ t^2/2 & \text{ if }\vert t\vert \le \lambda \end{cases} \]

3.8.3 Dantzig 选择器

Candes and Tao (2007)提出准则： \[ \text{min}_\beta\Vert\beta\Vert_1\text{ subject to }\Vert \mathbf X^T(\mathbf y-\mathbf X\beta)\Vert_\infty\le s\tag{3.78} \] 这个解是 Dantzig selector (DS). 等价式是 \[ \min_\beta \Vert \mathbf X^T(\mathbf y-\mathbf X\beta)\Vert_\infty\text{ subject to } \Vert\beta\Vert_1\le t\tag{3.79} \]

Candes and Tao (2007)证明了求解 DS 是线性规划问题 Dantzig 选择器可以最小化当前残差与所有变量之间的最大内积.

3.8.4 The Grouped Lasso

在一些预测变量属于预定义的群体的情形中；比如属于同一个生物路径的基因, 或者表示类别型数据层次的指示变量. 要对群体中每个成员一起进行收缩或选择. Grouped lasso 便可以实现.

假设 \(p\) 个预测变量被分到 \(L\) 个群中, 在第 \(\ell\) 个群中有 \(p_\ell\) 个成员. 采用矩阵 \(\mathbf X_\ell\) 来表示对应第 \(\ell\) 个群的预测变量, 对应的系数向量为\(\beta_\ell\). grouped-lasso 最小化下面的凸准则 \[ \underset{\beta\in \mathbb{R}^p}{\text{min}}\Big(\Vert \mathbf y-\beta_0\boldsymbol 1-\sum\limits_{\ell=1}^L\mathbf X_\ell\beta_\ell\Vert_2^2+\lambda \sum\limits_{\ell=1}^L\sqrt{p_\ell}\Vert \beta_\ell\Vert_2\Big)\tag{3.80} \] 在保持模型预测准确性的同时实现特征选择和正则化.提高模型可解释性. ## 3.8.5 lasso 的更多性质

当 \(N\) 和 \(p\) 增长时, 研究lasso 的能力以及重建模型相关的过程.例如, Donoho (2006b)集中在 \(p>N\) 的情形并且当边界 \(t\) 变大时 lasso 的解. 取极限得出在所有零训练误差的模型中最小的 \(L_1\) 范数解.

许多的相关结果对模型矩阵假设了如下条件 \[ \underset{j\in \mathcal S^c}{\text{max}}\Vert \mathbf x_j^T\mathbf X_{\mathcal S}(\mathbf X_{\mathcal S}^T\mathbf X_{\mathcal S})^{-1}\Vert_1\le (1-\epsilon)\text{ for some }\epsilon\in (0, 1]\tag{3.81} \]

修改 lasso 惩罚函数使得更大的系数收缩得不要太剧烈； 平稳削减绝对偏差法 用 \(J_a(\beta,\lambda)\) 替换 \(\lambda\vert\beta \vert\), 其中 \(a\ge 2\) \[ \frac{dJ_a(\beta,\lambda)}{d\beta}=\lambda \cdot \text{sign}(\beta)[I(\vert \beta\vert\le \lambda)+\frac{(a\lambda-\vert \beta\vert)_+}{(a-1)\lambda}I(\vert\beta\vert>\lambda)]\tag{3.82} \]

3.8.6 Pathwise Coordinate Optimization

简单坐标下降是一种替代计算 lasso 的 LARS 算法. 是固定 Lagrangian 形式( 3.52 ) 中的惩罚参数 \(\lambda\), 在控制其它参数固定不变时, 相继地优化每一个参数.

note “Recall”

\[ \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} \]

假设预测变量都经过标准化得到 0 均值和单位范数. 用 \(\tilde\beta_k(\lambda)\) 表示惩罚参数为 \(\lambda\) 时对 \(\beta_k\) 的当前估计. 我们可以分离出 式( 3.52 ) 的 \(\beta_j\),

\[ F(\tilde\beta(\lambda),\beta_j)=\frac{1}{2}\sum\limits_{i=1}^N(y_i-\sum\limits_{k\neq j}x_{ik}\tilde \beta_k(\lambda)-x_{ij}\beta_j)^2+\lambda \sum\limits_{k\neq j}\vert \tilde \beta_k(\lambda)\vert+\lambda \vert \beta_j\vert \]

压缩了截距并且为了方便引出因子 \(\frac 12\). 看成是响应变量为部分残差 \(y_i-\tilde y_i^{(j)}=y_i-\sum_{k\neq j}x_{ik}\tilde \beta_k(\lambda)\). 得到下面的等式 \[ \tilde \beta_j(\lambda)\leftarrow S(\sum_{i=1}^Nx_{ij}(y_i-\tilde y_i^{(j)}),\lambda)\tag{3.84} \]

\(S(t,\lambda)=\text{sign}(t)(\vert t\vert-\lambda)_+\) 是表 3.4 中的软阈值算子. \(S(\cdot)\) 中的第一个变量是部分残差在标准化变量 \(x_{ij}\) 上的简单最小二乘系数. 式( 3.84 ) 的重复迭代得到 lasso 估计 \(\hat\beta(\lambda)\).

这种算法可以有效地计算在 \(\lambda\) 的每个网格结点上 lasso 的解.

3.9 计算上的考虑

通过 LAR 算法实现的 lasso 的计算量与最小二乘拟合有相同的阶数.

“weiya 注：Cholesky and QR decomposition for Least Squares”

对于 Cholesky 分解, 考虑线性方程 \(X^TX\beta=X^TY\), 若 \(X^TX=LL^T\), 则转换为求解 \(L w=X^TY\) and \(L^T\beta=w\).

对于 QR 分解, 考虑线性方程 \(X \beta = Y\), 若 \(X = QR\), 则 \(QR \beta = Y=QQ^TY\), 所以转换为求解 \(R \beta=Q^TY\).