ESL读书笔记：Chap 3 Linear Methods for Regression

Chap 3 Linear Methods for Regression

Introduction

• 线性回归是假设$E(Y|X)$是关于X的线性组合
• 线性模型的可解释性强，在数据量少、信噪比低或稀疏的数据上往往可以得到比非线性模型更好的结果
• 通过对X进行变换，很容易达到非线性拟合的效果，这种衍生方法有时也被称作basis-function methods（基函数方法）

Linear Regression Models and Least Squres

从多个角度理解OLSE

• 在不加任何分布假设，只是假设$Y=f_{\theta }(X)+ϵ$下进行。在多元线性回归模型中，解为：
  $$\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$$
• 在高斯-马尔可夫假设下，OLSE是BLUE的，即对于任何$\alpha \in {\Re }^{p\times 1}$，在${\alpha }^T\beta$的所有关于Y的线性无偏估计中，${\alpha }^T\hat{\beta }$是方差最小的；$\hat{V}=Var(\hat{\beta })$,$\tilde{V}=Var(\tilde{\beta })$，其中$\tilde{\beta }$是$\beta$的任一无偏线性估计，则有$\hat{\mathbf{V}}⪯\tilde{\mathbf{V}}.$，即$\tilde{V}-\hat{V}$是半正定的。

【注】：在证明中使用$e$与$\hat{\beta }$独立即可

• $\hat{Y}=X\hat{\beta }$实际上是将Y投影到X组成的超平面，残差e与$\hat{Y}$正交

Gram-Schmidt procedure for multiple regression

该方法从另一个角度去理解多元线性回归模型的系数估计，并且提供了一种新的数值估计的方法。

首先需要明确的是，如果X输入变量之间是正交的，那么多元线性模型系数估计的值和单变量模型估计的值得到的值是一样的。

【注】：但是由于估计的方差是不一样的，所以检验是不一样的

由于正交的设计矩阵在计算上有很好的方便，因此一种很自然的想法是希望调整X使得是正交的，也就是接下来介绍的Gam-Schimidt算法。首先只看一维的情况来获得一个直观的认识。假如$X\in {\Re }^{N\times 2}=[1,X_1]$，那么估计：
$${\hat{\beta }}_1=\frac{⟨\mathbf{x}-\bar{x}1,\mathbf{y}⟩}{⟨\mathbf{x}-\bar{x}1,\mathbf{x}-\bar{x}1⟩},$$
其中$1$可以认为是$X_0$，上述过程可以认为是以下步骤：

• 使用$X_0$对$X_1$进行线性回归，得到残差$z=X_1-\bar{x}1$，易知残差均值为0
• 使用残差$z$对Y进行线性回归，系数估计值恰好为${\hat{\beta }}_1$

延续这个想法，我们可以得到下面的Gam-Schimidt算法:

从这个算法可以看出：

• $x_j$是$z_0,...,z_j$的线性组合

事实上我们有：
$$\begin{array}{rl}X=& Z\Gamma \\ =& Z{D}^{-1}D\Gamma \\ =& QR\end{array}$$
其中$\Gamma$是一个上三角阵，$D$是一个对角阵，且$D_{jj}=||z_j|{|}_2$。X=QR是X的一个QR分解，$Q\in {\Re }^{N\times (p+1)}$是正交的$Q^{T}Q=I$，R是上三角阵$R\in {\Re }^{(p+1)\times (p+1)}$，从而有$\hat{\beta }=R^{-1}{Q}^{T}Y,\hat{Y}=Q{Q}^{T}Y$

• 基于$Z$进行系数估计得到的$\hat{\beta }$项中，只有${\hat{\beta }}_p$是和$x_p$有关的

所以综合来说基于Z估计得到的$\widehat{{\beta }_p}$就是基于X估计得到的$\widehat{{\beta }_p}$。 Gram-Schmidt方法实际上就是在此基础上进行的。

【注】：这个算法中由于$z_j$是$x_j$与$z_0,...,z_{j-1}$做回归得到的残差，所以$z_j$与$z_0,...,z_{j-1}$是正交的，最终做回归的那个$Z$也会是一个列正交的矩阵，所以对整个Z进行系数估计就等价于分别进行，使用一元线性回归，有利于计算。

多变量多元线性回归

相较于原始的单变量多元线性回归，多变量多元线性回归只是在原始的基础上对Y进行了拓展：
$$\begin{array}{rl}& Y\in {\Re }^{N\times 1}\to Y\in {\Re }^{N\times K}\\ & ϵ\in {\Re }^{N\times 1}\to E\in {\Re }^{N\times K}\\ & Y=X\beta +ϵ\to Y=XB+E\end{array}$$

损失函数可以是原始的RSS：
$$\begin{array}{rl}RSS(B)=& \sum _{k=1}^K\sum _{i=1}^N(y_{ik}-f_k(x_i){)}^2\\ =& tr[(Y-XB{)}^T(Y-XB)]\end{array}$$
易知系数估计为：
$$\begin{array}{r}\hat{B}=(X^{T}X{)}^{-1}{X}^{T}Y\end{array}$$
但是对于$E$不是对角阵的情况，基于GM假设，最大化似然函数，使用加权的RSS是更常见的情况，也就是：
$$\begin{array}{r}RSS(B,\Sigma )=\sum _{i=1}^N(y_i-f(x_i){)}^T{\Sigma }^{-1}(y_i-f(x_i))\end{array}$$
但实际上，不论中间的矩阵是什么，只要是二次型的损失函数，系数估计的结果都不变。

# 证明：
记OLSE所得系数估计为$\hat{B}$，对其余估计$B,f(x_i)=B^T{x}_i$有：
$$\begin{array}{rl}RSS(B,\Sigma )=& \sum _{i=1}^N(y_i-f(x_i){)}^T{\Sigma }^{-1}(y_i-f(x_i))\\ =& \sum _{i=1}^N(y_i-B^T{x}_i+{\hat{B}}^T{x}_i-{\hat{B}}^T{x}_i{)}^T{\Sigma }^{-1}(y_i-B^T{x}_i+{\hat{B}}^T{x}_i-{\hat{B}}^T{x}_i)\\ =& \sum _{i=1}^N(y_i-{\hat{B}}^T{x}_i{)}^T{\Sigma }^{-1}(y_i-{\hat{B}}^T{x}_i)(\ast )\\ & +\sum _{i=1}^N({\hat{B}}^T{x}_i-B^T{x}_i{)}^T{\Sigma }^{-1}({\hat{B}}^T{x}_i-B^T{x}_i)(\ast \ast )\\ & +2\sum _{i=1}^N(y_i-{\hat{B}}^T{x}_i{)}^T{\Sigma }^{-1}({\hat{B}}^T{x}_i-B^T{x}_i)(\ast \ast \ast )\end{array}$$

其中$(\ast )=RSS(\hat{B},\Sigma ),(\ast \ast )\ge 0$，现在证明$(\ast \ast \ast )=0$.

若记${\alpha }_i=[0,...,1,\mathrm{0...},0{]}^T$则:
$$\begin{array}{r}{y}_i^T={\alpha }_i^{T}Y\\ x_i^T={\alpha }_i^{T}X\end{array}$$

从而
$$\begin{array}{rl}(\ast \ast \ast )=& \sum _{i=1}^N(y_i-{\hat{B}}^T{x}_i{)}^T{\Sigma }^{-1}({\hat{B}}^T{x}_i-B^T{x}_i)\\ & =\sum _{i=1}^N{\alpha }_i^T(I_N-H)Y{\Sigma }^{-1}(\hat{B}-B{)}^T{X}^T{\alpha }_i\\ & :=\sum _{i=1}^N{\alpha }_i\Gamma {\alpha }_i^T\\ & =trace(\Gamma )\\ & =trace((I_N-H)Y{\Sigma }^{-1}(\hat{B}-B{)}^T{X}^T)\\ & =trace(X^T(I_N-H)Y{\Sigma }^{-1}(\hat{B}-B{)}^T)\\ & =0\end{array}$$
从而原命题成立.#

补充

• 对${\hat{\beta }}_i$进行区间估计的时候，就算高斯-马尔可夫假设不成立，在$N\to +\infty$时原始区间估计仍然有近似区间成立
• 在GM假设下，预测新点$x_0$有：
• $$\begin{array}{rl}E(Y_0-\hat{f}(x_0){)}^2=& {\sigma }^2+E(x_0^T\hat{\beta }-f(x_0){)}^2\\ =& {\sigma }^2+MSE(\hat{f}(x_0))\end{array}$$

Subset Selection

全模型OLSE存在的问题：

• 虽然是无偏的，但是方差大；
• 过多变量不利于解释

因此需要使用选模型，方法有forward-stepwise selection, backward-step wise selection,forward- and backward-stepwise selection。这里补充一个新的方法：forward-stagewise selection。这个方法相对来说低效，但是在高维的时候有很好的效果。

Forward-stagewise regression (FS)：

• 初始化：X中心化，使用$1$对Y进行回归得到$\hat{Y}=\bar{Y}$，得到残差$R_0=Y-\bar{Y}$
• 循环迭代：第k步都有基于前k-1步得到的残差$R_{k-1}$，找到X中与$R_{k-1}$最相关的变量进行单变量回归系数${\hat{\beta }}_k$，并在不改变先前变量的系数估计条件下加入模型，得到残差$R_k$
• 直到没有变量与当前残差相关（相关性高于某个阈值）停止迭代

与stepwise方法最不相同之处在于stagewise每引入一个新的变量，不会改变之前引入变量的系数估计，针对残差进行一元线性回归。也因此，stagewise可能会进行大于p步的迭代，这样的slow fitting在解决高维问题很有优势。

Shrinkage Methods

变量子集选择相对来说是一个离散的过程，自变量要么纳入回归要么完全不考虑，这可能会导致高方差的问题，而Shrinkage Methods则是用一种连续的方式来对变量进行“选择”。

Ridge Regression

$${\hat{\beta }}^{ridge}=\underset{\beta }{\mathrm{argmin}}\left\{\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^p{x}_{ij}{\beta }_j{)}^2+\lambda \sum _{j=1}^p{\beta }_j^2\right\}$$

【注】：在神经网络中引入这样的惩罚项也称作weight decay

以上也等价于：
$$\begin{array}{rl}& {\hat{\beta }}^{ridge}=\underset{\beta }{\mathrm{argmin}}\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^p{x}_{ij}{\beta }_j)\\ & s.t.\sum _{j=1}^p{\beta }_j^2\le t\end{array}$$

需要注意的是：

• 惩罚项中并没有对${\beta }_0$进行处理！最后对${\beta }_0$的估计只需要${\hat{\beta }}_0={\sum }_{i=1}^N(y_i-{\sum }_{j=1}^p{x}_{ij}{\beta }_j)$
• 需要先进行数据的标准化！

求解得到估计为：

$${\hat{\beta }}^{ridge}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}Y$$

由此可见：

• 使用岭回归可以有效缓解多重共线性（$X^{T}X+\lambda I$是可逆的）
• 这等价于系数$\beta$先验分布为$N(0,{\tau }^2)$并且是相互独立的，在此情况下的后验均值，所以实际上和贝叶斯统计得到的解是一样的

接下来从另一个角度来理解岭回归：SVD分解+PCA

首先来介绍一下SVD分解：

$$X\in {\Re }^{N\times p}:X=UD{V}^T$$

其中，$U\in {\Re }^{N\times p},V\in {\Re }^{p\times p}$是正交的，D是$p\times p$的对角阵，并且对角元素递减：$d_1\ge d_2\ge ...\ge d_p\ge 0$

那么此时用SVD分解代替表示X有：
$$X{\hat{\beta }}^{OLS}=U{U}^{T}Y$$
$$X{\hat{\beta }}^{ridge}=UD(D^2+\lambda I{)}^{-1}D{U}^{T}Y=\sum _{j=1}^p{u}_j\frac{d_j^2}{d_j^2+\lambda }{u}_j^{T}Y$$
可以看到，估计$\hat{Y}$实际上都是Y的投影，但是相对来说，岭回归的估计值在原始的${\sum }_{j=1}{u}_j{u}_j^{T}Y$的基础上还进行了一个加权$\frac{d_j^2}{d_j^2+\lambda },\lambda \le 0$，也就是说会相对来说缩减为原始的$\frac{d_j^2}{d_j^2+\lambda }$。接下来结合PCA对这个缩减系数进行解释。事实上，对X进行SVD分解$X=UD{V}^T$，假设X已经中心化，那么主成分Z为$UD$，所以:
$${\hat{Y}}^{ridge}=Z(D^2\lambda I{)}^{(-1)}=\sum _{j=1}^p\frac{1}{d_j^2+\lambda }{z}_j{z}_j^{T}Y$$
相当于是在Z（未标准化标准差）的方向上投影再加权。

在对X进行标准化之后，实际上有：
$$S=X^{T}X/N=V{D}^2{V}^T/N$$
又注意到V是正交的，$D^2$是对角阵并且元素是递减的，因此可以认为$D^2/N$就是PCA中那个对角阵$\Lambda$！$z_j=X{v}_j=u_j{d}_j$，并且对于标准化后的X$var(z_j)=d_j^2/N$，因此$u_j$可以看成是标准化后的$z_j$。因此$u_j{u}_j^{T}Y$可以看成是在第j个主成分上的投影，$d_j^2$越大实际上说明对应的主成分承载的信息越多，方差越大，那么$\frac{d_j^2}{d_j^2+\lambda }$越大，所以可以解释为在承载信息更多的主成分方向上投影的越多。

进一步的，相对于OLSE可以定义岭回归的自变量自由度：
$$\begin{array}{rl}df(\lambda )=& tr(X(X^{T}X+\lambda I{)}^{-1}{X}^T)\\ =& tr(H_{\lambda })\\ =& \sum _{j=1}^p\frac{d_j^2}{d_j^2+\lambda }\end{array}$$

易知，$\lambda \to 0,df(\lambda )\to p;\lambda \to \infty ,df(\lambda )\to 0$。

岭回归的结果也可以理解为对X和Y分别进行增广。
$$\tilde{X}=\left[\begin{array}{c}X\\ \sqrt{\lambda \alpha }{I}_P\end{array}\right],\tilde{Y}=\left[\begin{array}{c}Y\\ 0_P\end{array}\right]$$
并使用OLSE得到$\hat{\beta }$不难验证就是岭回归得到的系数估计值，从直观上也不难理解会将$\beta$向0压缩。elastic net也可以使用增广去理解，相当于增广+Lasso。

Lasso

$${\hat{\beta }}^{lasso}=\underset{\beta }{\mathrm{argmin}}\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^p{x}_{ij}{\beta }_j{)}^2+\lambda \sum _{j=1}^p|{\beta }_j|$$

也等价于：
$$\begin{array}{rl}& {\hat{\beta }}^{lasso}=\underset{\beta }{\mathrm{argmin}}\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^p{x}_{ij}{\beta }_j{)}^2\\ & s.t.\sum _{j=1}^p|{\beta }_j|\ge t\end{array}$$
相较于岭回归，lasso不同的唯一点在于惩罚项由L2变成了L1。这在计算上并不利于求解，因为在0不可导；但是在解释上更好理解，当(30)式中的t等于OLSE对应的$||{\hat{\beta }}^{OLS}|{|}_1$则lasso的解和OLSE相同，而如果将t缩减为原始的1/2，则可以近似认为系数估计的值平均来说缩减为原始的1/2。

对比子集选择、岭回归、Lasso

对比子集选择、岭回归、lasso，三者对于变量选择都有一定的作用：

• 如果X标准化并且是正交的，那么以上三种方法都有显式解：

Estimator	Formula
Best subset(size M)	$\hat\beta_j I(
Ridge	${\hat{\beta }}_j/(1+\lambda )$
Lasso	$sign(\hat \beta_j)(

• 相对来说，Lasso可以真正的将某一个变量的系数估计变成0，做到变量选择，但是岭回归不能
• 将惩罚项进行拓展，$penaty(\lambda ,q)=\lambda {\sum }_{j=1}^p|{\beta }_j{|}^q$，当q小于1的时候，会在轴向上有更多的空间：

• 在p=1的时候，0是不可导的，这个特性比较特殊，所以引入一种新的正则方法：elastic penalty，对应的penalty为$\lambda L_2+(1-\lambda )L_1$

Least Angle Regression

这个网站写的比较清楚：https://keson96.github.io/2016/10/26/2016-10-26-Least-Angle-Regression/

LAR算法的过程大致如下（X,Y已经标准化）：

每一轮迭代的时候做的事情就是沿着多元OLS系数估计的角平分线方向进行小步幅的移动，这样会得到在相应移动之后的残差，由于并没有完全到达OLSE，所以此时的残差和“回归”变量之间并不是不相关的。一直移动直到开始出现另一个自变量与该残差的相关性与当前所有自变量与该残差的相关程度$|\rho |$相同，并开始进行下一轮的新迭代。

【注】：从直观的角度去理解，shrinkage methods都是基于原始估计值在“缩减”，LAR就是不停地在OLSE方向上移动小步幅，直到所有active set中的自变量和残差的相关程度相同。

具体来说，在每一轮更新系数估计的时候，在第k轮，${\mathcal{A}}_k$为当前的k-1个变量，${\beta }_{{\mathcal{A}}_k}$为当前还没有引入第k个变量时的系数估计，残差$r_k=y-X_{{\mathcal{A}}_k}{\beta }_{{\mathcal{A}}_k}$，那么下一步的更新方向即为：
$${\delta }_k=(X_{{\mathcal{A}}_k}^T{X}_{{\mathcal{A}}_k}{)}^{-1}{X}_{{\mathcal{A}}_k}^T{r}_k$$
更新系数为：
$${\beta }_{{\mathcal{A}}_k}(\alpha )={\beta }_{{\mathcal{A}}_k}+\alpha {\delta }_k$$

其中$\alpha$是一个很小的数。一直调整系数估计${\beta }_{{\mathcal{A}}_k}$至有一个新的变量和当前变量的回归残差的相关程度达到某一个阈值要求，就引入这个新的变量，并重复上述更新过程。所以实际上在第k步，如果实现已经进行了标准化，用inner product计算相关系数有：
$$x_j^T(Y-X\beta )=\gamma s_j,\forall x_j\in {\mathcal{A}}_k$$

其中$\gamma$是相关系数的绝对值，$s_i=sign(x_i^T\ast (Y-X\beta ))$，也就是说当前${\mathcal{A}}_k$中所有变量和当前残差的相关程度是相同的。具体的这幅图可以很好地反映。

我们可以对比随着LAR迭代$||\hat{\beta }|{|}_1$增加，各系数的变化路径，和Lasso随着$\lambda$减小各系数的变化路径，可以看到二者唯一不同的地方在于某一个系数变为0之后，Lasso会直接将其踢出当前回归器，而LAR会接着朝某一方向移动。因此参照Lasso进行调整。

**LAR和Lasso结合的方法：当系数调整时等于0，则舍去改变量，提出当前active set ${\mathcal{A}}_k$！**这个方法和Lasso的结果相近，并且计算很高效，尤其是$p>>N$；只需要p步就行，而lasso往往需要大于p步。

**LAR和Lasso为什么在效果上这么相似？**我们从数学推导更严谨的角度上，直接关注最终active set的结果。假设都已经标准化，使用inner-product来度量相关程度。

对于LAR，假设最终的active set是$\mathcal{A}$，那么有：
$$x_j^T(Y-X\beta )=\gamma s_j,\forall x_j\in \mathcal{A},s_j=sign(x_j^T(Y-X\beta ))$$
$$|x_j^T(Y-X\beta )|\le \gamma ,\forall x_j\notin \mathcal{A}$$
其中$\gamma$就是那个共同的相关程度。

而对于Lasso，如果是真正的active set$\mathcal{B}$，有${\beta }_j\ne 0$，所以对于active set中的变量进行系数估计：
$$R(\beta )=1/2||Y-X\beta |{|}_2^2+\lambda ||\beta |{|}_1$$

实际上是可导的，并且导数等于零有：

$$x_j^T(Y-X\beta )=\lambda sign({\beta }_j),\forall j\in \mathcal{B}$$
而对于non-active变量自然会有$x_j^T(Y-X\beta )\le \lambda ,\forall j\notin \mathcal{B}$。

自由度的定义：LAR & Lasso

重新定义自由度：对于N个观测$Y\in {\Re }^{N\times 1}$，根据模型得到估计$\hat{Y}=({\hat{y}}_1,...,{\hat{y}}_N)$，并定义模型自由度为，
$$df(\hat{y})=\frac{1}{{\sigma }^2}\sum _{i=1}^{N}Cov({\hat{y}}_i,y_i)$$

1. 容易验证，在此定义下，全模型的自由度为p+1；岭回归的自由度为$tr(H_{\lambda })$
2. 规定选模型选择k个自变量，使用随机模拟容易验证df会大于k
3. 对于LAR，第k步的时候，自由度为k
4. 对lasso+LAR，自由度近似为最终active set的变量个数

Methods Using Derived Input Directions

在很多时候可能会出现多重共线性，有一种解决方法是对原始输入$X_j$做线性组合提取$z_1,...,z_M$，并使用$Z$代替X进行回归，最后再变换回X。

PCR

PCR：主成分分析+回归。与岭回归进行对比，从之前的铺垫（X的SVD分解）可以知道，岭回归是一种shrinkage methods，对于方差较小的主成分，会做一些缩减，而对比来说，PCR则是对于方差很小的主成分直接缩减到0，较大的主成分则不做任何缩减，如图。

偏最小二乘回归

偏最小二乘回归和PCR相同点在于都是对X进行线性组合得到主成分Z，再使用Z进行回归，最后还原回使用X的回归。\textbf{不同点在于PCR得到主成分实际上只考虑了X，而偏最小二乘回归还考虑到了主成分Z和Y之间的相关性}，并且适用于多对多的建模问题。当两组变量的个数很多，并且都存在多重共线性，而观测数据的数量又较少的时候可以使用。

现在考虑$Y_2,...Y_p$与自变量$X_1,...,X_m$的建模问题。偏最小二乘回归的基本做法是：

• 首先在自变量集中提取第一成分$T_1$，$T_1$是X的线性组合，并且能够尽可能多地提取X的变异信息；

• 同时在因变量集Y中提取第一成分$U_1$，并且要求$T_1$与$U_1$相关程度达到最大，建立因变量与$T_1$的回归，如果回归的精度已经达到要求则停止，否则继续；

• 如果精度不符合要求，则构建$X$与$T_1,...,T_k$的回归方程，得到残差$E_k$，将$E_k$作为自变量集；构建$Y$与$T_1,...,T_k$的回归方程，得到残差$F_k$，将$F_k$作为因变量集，提取第k+1个成分$T_{k+1}$；
  $$t_1=X{w}_1$$
  $$u_1=Y{v}_1$$
  假设X和Y都已经标准化，现要求同时最大化$T_1$的方差，$U_1$的方差和$T_1$与$U_1$的相关程度。
  $$\underset{w_1,v_1}{\max }\sqrt{Var(T_1)Var(U_1)}\rho (T_1,U_1)=\underset{w_1,v_1}{\max }Cov(T_1,U_1)=w_1^T{X}^{T}Y{v}_1$$
  使得：
  $$||w_1|{|}_2=1,||v_1|{|}^2=1$$
  解得，$w_1$是$X_0^T{Y}_0{Y}_0^T{X}_0$最大特征根对应的特征向量，$v_1$是$Y_0^T{X}_0{X}_0^T{Y}_0$最大特征根对应的特征向量。

• 若最终得到对自变量集提取的r个主成分$T_1,...,T_r$，偏最小二乘回归建立$T_1,...,T_r$对$Y_1,...,Y_p$的回归，再重新表示为原始自变量与Y的回归即可

Discussion:A Comparison of Selection and Shrinkage Methods

• 岭回归效果最好的那一个。它对于所有主成分的方向都会做缩减，但是尤其在低方差的方向缩减，并且是一个相对连续光滑的shrinkage
• 偏最小二乘回归在低方差方向做缩减，同时在高方差方向可能做膨胀，这使得PLS有一点不稳定

实际上就是，PLS、子集选择、PCA都是离散化的shrinkage过程，而岭回归、Lasso则是连续的。

Multiple Outcome Shrinkage and Selection

对于有多个输出的$Y\in {\Re }^{N\times K}$中，如何使用shrinkage。一种想法是分开$Y_1,Y_2,...,Y_K$，分别使用相同超参数或不同超参数下的Lasso或岭回归。另一种想法则是延续PCA方法的canonical correlation analysis(CCA)，本节介绍CCA。

CCA希望能够找到一组X的线性组合$X{v}_m,m=1,\mathrm{2...}M$和相应的Y的线性组合$Y{u}_m,m=1,\mathrm{2...}M$使得：
$${\mathrm{Corr}}^2(\mathbf{Y}{u}_m,\mathbf{X}{v}_m)$$
最大。注意M至多取min(K,p).

CCA的解：假设X,Y均已中心化

最大化：
$${\mathrm{Corr}}^2(\mathbf{Y}{u}_m,\mathbf{X}{v}_m)$$
等价于:
$$\underset{\begin{array}{c}{u}^T({\mathbf{Y}}^T\mathbf{Y})u=1\\ v^T({\mathbf{X}}^T\mathbf{X})v=1\end{array}}{\max }{u}^T({\mathbf{Y}}^T\mathbf{X})v,$$
如果先对$Y$做“标准化”，使得变量之间是不相关且等方差的，即$\tilde{Y}=Y(Y^{T}Y{)}^{-1/2}$，对X同理得到$\tilde{X}=X(X^{T}X{)}^{(-1/2)}$；并令$\tilde{u}=((Y^{T}Y{)}^{1/2})u,\tilde{v}=((Y^{T}Y{)}^{1/2})v$。那么上述问题等价于：
$$\underset{\begin{array}{r}{\tilde{u}}^T\tilde{u}=1\\ {\tilde{v}}^T\tilde{v}=1\end{array}}{\max }{\tilde{u}}^T(\widetilde{{\mathbf{Y}}^T}\tilde{\mathbf{X}})\tilde{v},$$
使用拉格朗日乘子法，可以得到：
$$\begin{array}{rl}M{\tilde{v}}_1& =\lambda {\tilde{u}}_1,\\ M^T{\tilde{u}}_1& =\lambda {\tilde{v}}_1.\end{array}$$
其中$M={\tilde{Y}}^T\tilde{X}$。对M使用SVD分解：
$$M=({\mathbf{Y}}^T\mathbf{Y}{)}^{-\frac{1}{2}}({\mathbf{Y}}^T\mathbf{X})({\mathbf{X}}^T\mathbf{X}{)}^{-\frac{1}{2}}=\tilde{U}D{\tilde{V}}^T$$
\textbf{使用代数相关知识}，可以得到${\tilde{u}}_1$是$\tilde{U}$的第一列，${\tilde{v}}_1$是$\tilde{V}$的第一列。那么$u=\tilde{u}((Y^{T}Y{)}^{-1/2}),v=\tilde{v}((X^{T}X{)}^{-1/2})$.

后续的成分提取，和PCA一致，需要与前面的成分不相关，即：
$$\underset{\begin{array}{c}{u}^T({\mathbf{Y}}^T\mathbf{Y})u=1\\ v^T({\mathbf{X}}^T\mathbf{X})v=1\\ u^T{u}_j=0,\forall j<k\\ v^T{v}_j=0,\forall j<k\end{array}}{\max }{u}^T({\mathbf{Y}}^T\mathbf{X})v,$$
不难证明得到：\textbf{对$M=({\mathbf{Y}}^T\mathbf{Y}{)}^{-\frac{1}{2}}({\mathbf{Y}}^T\mathbf{X})({\mathbf{X}}^T\mathbf{X}{)}^{-\frac{1}{2}}=\tilde{U}D{\tilde{V}}^T$进行SVD分解，$u=\tilde{u}((Y^{T}Y{)}^{-1/2}),v=\tilde{v}((X^{T}X{)}^{-1/2})$，即可得到相应线性变化的解。}

此外，对于Reduced-rank regression，假设$Var(ϵ)=\Sigma$，那么希望解决：
$${\hat{\mathbf{B}}}^{\text{rr}}(m)=\underset{\text{rank(B)}=m}{\mathrm{argmin}}\sum _{i=1}^N(y_i-{\mathbf{B}}^T{x}_i{)}^T{\mathbf{\Sigma }}^{-1}(y_i-{\mathbf{B}}^T{x}_i).$$
如果$\Sigma =Y^{T}Y/N$,那么可以使用CCA给出解：
$${\hat{\mathbf{B}}}^{\mathrm{rr}}(m)=\hat{\mathbf{B}}{\mathbf{U}}_m{\mathbf{U}}_m^{-},$$
其中$\hat{B}$是最小二乘下得到的系数估计，$U_m=[u_1,...u_m]$，$U_m^{-}={\Sigma }^{-1/2}{\tilde{U}}_m$。由此：
$${\hat{\mathbf{B}}}^{\mathrm{rr}}(M)=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T(\mathbf{Y}{\mathbf{U}}_m){\mathbf{U}}_m,$$
可以理解为是使用X对$Y{U}_m$进行回归，再使用$U_m^{-}$回到原本的系数上。

预测：
$$\begin{array}{rl}{\hat{\mathbf{Y}}}^{\mathrm{rr}}(m)& =\mathbf{X}({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{Y}{\mathbf{U}}_m{\mathbf{U}}_m^{-}\\ & =\mathbf{HY}{\mathbf{P}}_m,\end{array}$$
其中H是OLSE的帽子矩阵，$P_m$是基于CCA得到的秩为m的矩阵。此外由于在证明过程中并没有对$\Sigma$做限制，所以以上形式是通用的。

基于reduced rank regression的一种shrinkage为：
$$\widehat{{\mathbf{B}}^{\mathbf{c}+\mathbf{w}}}=\widehat{\mathbf{BU\Lambda }{\mathbf{U}}^{-1}},$$
其中$\Lambda$是一个对角阵，可以使用CV方法得到。由此，预测为：
$$\widehat{{\mathrm{Y}}^{\mathrm{c}+\mathrm{w}}}=\mathrm{HY}{\mathrm{S}}^{\mathrm{c}+\mathrm{w}},$$
更一般的，加上岭回归：
$$\widehat{{\mathbf{Y}}^{\mathrm{ridge},\mathrm{c}+\mathrm{w}}}={\mathbf{A}}_{\lambda }\mathbf{Y}{\mathbf{S}}^{\mathbf{c}+\mathbf{w}},$$
$${\mathbf{A}}_{\lambda }=\mathbf{X}({\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^T$$

More on the Lasso and Related Path Algorithms

Incremental Forward Stagewise Regression

首先回顾一下Lasso和LAR。从前面的推导我们可以发现，对于最终的active set中的变量，Lasso的解需要满足：
$$x_j^T(Y-X\beta )=\lambda sign({\beta }_j),\forall j\in \mathcal{B}$$
而对于non-active变量自然会有$x_j^T(Y-X\beta )\le \lambda ,\forall j\notin \mathcal{B}$。所以就会发现，此时残差与active set中的变量的相关程度是一致的！ 这也和LAR的想法是一样的，所以才可以使用LAR作为Lasso的一个近似求解。

现在想讨论的是，使用类似这种迭代算法（例如LAR）进行求解，达到shrinkage的目的，还有什么可以改进的，以及使用这种迭代算法，迭代的路径是什么。

本节首先介绍一个新的方法：incremental forward stagewise regression——$F{S}_{ϵ}$，并且和之前提到的stagewise regression进行对比

我们将它和LAR以及stagewise regression进行对比：

与LAR对比，直接关注到迭代过程，可以看到：

• LAR实际上每一次迭代是对当前active set中的所有系数估计进行迭代，迭代的方向是当前残差和active set中自变量的系数（针对多元线性模型）估计值方向，以使得残差和自变量越来越不相关，直到non-active set中出现和这些变量相关程度一样的新变量引入active set；

• 而FS算法则是不会更新active set中已有的变量，而是平等地对所有变量进行考察，对与当前残差最相关的变量$x_j$移动，移动的方向与一元线性模型估计得到的系数估计方向一致（这就和LAR不一样了）

• 与stagewise regression对比，一般的stagewise regression与$F{S}_{ϵ}$唯一不同的在于，更新的步幅并不是实现指定好的$ϵ$，而是$<x_j,r>$

简单来说，LAR每一次是缓慢移动所有的变量系数估计，而FS则是只动一个；stagewise每一次把${\beta }_i$移动拉满，而FS则只是移动一小步。

从LAR与FS的对比可以看到，LAR在迭代系数的时候，由于考虑的是一个多元线性模型的系数估计方向，所以可能和单独的一元线性回归的系数估计方向相反，于是可以做出以下的调整，以强制使得和一元线性回归估计的系数方向一致：

LAR、lasso、$F{S}_{0}modification$都可以使用LAR进行迭代计算，他们的路径都是分段线性的：

• 对于$ϵ\to 0$的情况$F{S}_0$，它和Lasso不同但都是分段线性的路径，书中提到它实际上也是在做一个优化："$F{S}_0$ is optimal per unit increase in L1 arc-length traveled along the coefficient path.Hence its coefficient path is discouraged from changing directions too often."在$p>>N$的时候相对于lasso，使用$F{S}_0$会更加smooth，效果更好;
• 当LAR每次迭代的时候都是单调的，三者都将是一致的；当不是单调的，但是不跨过0，LAR和lasso是一致的

Piecewise-Linear Path Algorithms

对于正则化的问题，都可以使用类似LAR的算法进行求解，都会有一个迭代路径，这样的算法称作"path algorithms"。现在关注这些正则条件对应的迭代路径。

$$\hat{\beta }(\lambda )=argmi{n}_{\beta }[R(\beta )+\lambda J(\beta )]$$
$$R(\beta )=\sum _{i=1}^{N}L(y_i,{\beta }_0+\sum _{j=1}^p{x}_j{\beta }_j)$$
其中L和J都是凸函数，则当：

• R是二次或分段二次的
• J是分段线性的

则以上迭代路径是分段线性的。

所以显然，$L_1,L_{\infty }$作为penalty在MSE下的迭代路径都是分段线性的。

The Dantzig Selector

$$\underset{\beta }{\min }||X^T(Y-X\beta )|{|}_{\infty }s.t.||\beta |{|}_1\le t$$
可以看到他是想让MSE的各方向梯度的最大值充分小，同时加上L1正则项。DS最终求解转化为一个线性规划问题。

不难理解这很像lasso：

• 在p<N，t充分大的时候，二者的解都会是OSLE；在p>N的时候，他们都会在L1正则项下得到一个解，而如果t较小，DS的解路径和Lasso不一样;
• Lasso是想让残差和active set中的所有变量相关程度相同，并且大于non-active set；DS则是想让最大的那一个相关程度最小，控制的是upper bound，但是这就有可能导致active set中的变量的相关程度有可能小于non-active set中的变量，这显然是不合的，所以DS有的时候并不准确

The Grouped Lasso

有的时候需要对一批变量同时进行正则，例如虚拟变量，于是出现the grouped lasso。
$$\underset{\beta \in {\Re }^p}{\min }(||Y-\beta 1-\sum _{l=1}^L{X}_l{\beta }_l|{|}_2+\lambda \sum _{l=1}^L\sqrt{p_l}||\beta )l|{|}_2)$$
其中$\sqrt{p_l}$表示第l组变量的变量个数。

Further Properties of the Lasso

对于Lasso来说，active set中的变量系数也会被shrink，一般来说并不是相合估计。一种减轻这种情况的方法是，在针对所有的变量使用Lasso筛选得到非零系数变量集之后，再针对Y使用一般的线性回归不用Lasso；另一种方法则是，在筛选得到非零系数变量集之后，重新使用CV得到${\lambda }^{\prime }<\lambda$，由此正则化的程度减轻，这种方法更加常用，被称作relaxed lasso.

另一种方法则是直接改变正则化的函数：
$$\frac{d{J}_a(\beta ,\lambda )}{d\beta }=\lambda \cdot \mathrm{sign}(\beta )[I(|\beta |\le \lambda )+\frac{(a\lambda -|\beta |{)}_{+}}{(a-1)\lambda }I(|\beta |>\lambda )]$$
相对Lasso来说，这个正则函数会对$|\beta |$较大的时候惩罚力度轻一些，但是并不是凸函数，计算上不方便。一种改进的想法（adaptive lasso）则是使用：
$$\sum _{j=1}^{\mathbf{p}}{w}_j|{\beta }_j|\mathrm{where}{w}_j=1/|{\hat{\beta }}_j{|}^{\nu },$$
其中${\hat{\beta }}_j$是使用OLSE得到的，$\nu >0$。adaptive lasso得到的系数估计具有相合性，并且保留凸性。

Pathwise Coordinate Optimization

另一种得到Lasso解的方法是simple coordinate descent.该方法的思想就是，每一次迭代的时候，轮流固定其他系数估计不变，针对某一个变量更新。假设所有X都标准化了，那么：
$$R(\tilde{\beta }(\lambda ),{\beta }_j)=\frac{1}{2}\sum _{i=1}^N{\left(y_i-\sum _{k\ne j}{x}_{ik}{\tilde{\beta }}_k(\lambda )-x_{ij}{\beta }_j\right)}^2+\lambda \sum _{k\ne j}|{\tilde{\beta }}_k(\lambda )|+\lambda |{\beta }_j|,$$
更新：
$${\tilde{\beta }}_j(\lambda )\leftarrow S(\sum _{i=1}^N{x}_{ij}(y_i-{\tilde{y}}_i^{(j)}),\lambda ).$$
其中$S(t,\lambda )=\text{sign}(t)(|t|-\lambda {)}_{+}$.此处S的前一项恰好是$x_j$对$Y-{\tilde{Y}}^{(j)}$进行回归得到的系数。对$j=1,\mathrm{2...}p$重复上述迭代过程，直到收敛，从而得到$\hat{\beta }(\lambda )$。

同样可以使用以上算法针对不同的进行。从最大的${\lambda }_{max}$，$\hat{\beta }({\lambda }_{max})=0$，每一次都减小一点点得到${\lambda }_k$，使用上述算法得到$\hat{\beta }({\lambda }_k)$，并使用上一个得到的$\hat{\beta }({\lambda }_{k-1})$作为初始值。这会比LARS更快，因为每一轮都是进行一元线性回归，计算会快很多。