7.3 正则化与高维线性回归

回顾我们上一章提到的最小二乘, 其中一个重要的假设即为设计矩阵 $\vec X$ 满列秩. 即样本数量 $n$ 和特征数量 $m$ 应满足 $n\geq m+1$ 的关系. 但是在很多数据中, 往往有着数以万计的特征却缺乏足够多的样本. 在此情况下 $n < m+1$ , 那么此时就有一些问题出现了. 比如此时方程的数量远远小于未知数的数量, 我们知道此时的方程便有无穷多组解. 此时的回归问题被称作高维线性回归 (High-dimensional Regression). 此时为了能够求解, 我们需要通过正则化 (Regularization) 来对我们的最小二乘进行适当的改造.

*{Ridge回归}

我们首先提出的是Ridge 回归 (Ridge Regression). Ridge回归在最小二乘的基础上加入了一个惩罚项 (Penalty Term), 使得改造过的矩阵重新变得可逆. 我们设 $\lambda>0$ 为一常数, 那么Ridge回归中参数向量 $\hat{\vec\beta}_{R,\lambda}$ 的解为

\hat{\vec\beta}_{R,\lambda} := \arg\min_{\vec\beta \in \mathbb{R}^{m+1}}\Big( ||\vec y - \vec X\vec\beta||^2+\lambda||\vec\beta||^2\Big).\tag{7.3}

推论 7.2

对任意的 $\vec X \in M_{mn}(\mathbb{R})$ 和参数 $\lambda>0$ , 矩阵 $\vec X^\top\vec X+\lambda\vec I$ 可逆.

证明

回顾第四章第 $4$ 节的知识, 对任意的 $\vec X$ 而言, $\vec X^\top\vec X$ 为半正定矩阵, 即 $\forall \vec u \in \mathbb{R}^n, \vec u^\top (\vec X^\top\vec X)\vec u \geq 0$ . 此时

\vec u^\top(\vec X^\top\vec X + \lambda\vec I)\vec u = \vec u^\top\vec X^\top\vec X\vec u + \lambda\vec u^\top\vec u>0,

因此 $\vec X^\top\vec X + \lambda\vec I$ 为正定矩阵, 而正定矩阵必然可逆.

∎

在Ridge回归下, 尽管设计矩阵 $\vec X$ 不为列满秩, 但 $\vec X^\top\vec X + \lambda\vec I$ 便为可逆矩阵, 我们也就可以在此情况下求出Ridge回归下的参数向量 $\hat{\vec\beta}_{R,\lambda}$ .

定理 7.6

若 $\lambda>0$ , 那么在(7.3)中, Ridge回归的参数向量有唯一解:

\hat{\vec\beta}_{R,\lambda} = \Big(\vec X^\top\vec X+\lambda\vec I \Big)^{-1}\vec X^\top\vec y.

证明

由定义可知,

\begin{aligned} ||\vec y - \vec X\vec\beta||^2+\lambda||\vec\beta||^2 &= (\vec y - \vec X\vec\beta)^\top(\vec y - \vec X\vec\beta)+\lambda\vec\beta^\top\vec\beta\\ &=\vec y^\top\vec y -\vec y^\top\vec X\vec\beta - (\vec X\vec\beta)^\top\vec y+ (\vec X\vec\beta)^\top(\vec X\vec\beta)+\lambda\vec\beta^\top\vec\beta\\ &=\vec y^\top\vec y - 2\vec y^\top\vec X\vec\beta+\vec\beta^\top\vec X^\top\vec X\vec\beta+\lambda\vec\beta^\top\vec\beta. \end{aligned}

设 $F(\vec\beta) = \vec y^\top\vec y - 2\vec y^\top\vec X\vec\beta+\vec\beta^\top\vec X^\top\vec X\vec\beta+\lambda\vec\beta^\top\vec\beta$ . 其梯度为

\nabla F(\vec\beta) = -2\vec X^\top\vec y + 2\vec X^\top\vec X\vec\beta+2\lambda\vec\beta

利用严格凸函数的性质, $\nabla F(\vec\beta)=\vec 0$ 的点对应唯一的极小值, 因此我们有 $\vec X^\top\vec y = \vec X^\top\vec X\vec\beta + \lambda\vec\beta$ . 由于 $\vec X^\top\vec X+\lambda\vec I$ 可逆. 则

\hat{\vec\beta}_{R,\lambda} = \Big(\vec X^\top \vec X+\lambda\vec I\Big)^{-1}\vec X^\top\vec y.

∎

我们知道, 最小二乘的几何意义可以看作是响应向量 $\vec y$ 在设计矩阵 $\vec X$ 的列空间上的正交投影. 那么Ridge回归有没有什么几何意义呢? 根据Lagrange乘数法, (7.3) 中的公式与下面的完全等价:

\min_{\vec\beta} || \vec y - \vec X\vec\beta||^2, \quad ||\beta||^2 \leq \lambda^2.

即Ridge回归便可看作是在普通最小二乘的前提下加上了关于参数向量 $\vec\beta$ 的约束条件. 为了方便理解, 我们先考虑只有两个参数的情况: 即 $\vec\beta = \begin{pmatrix} \beta_1 & \beta_2 \end{pmatrix}^\top$ . 此时

||\vec y - \vec X\vec\beta||^2 = (\vec y-\vec X\beta)^\top(\vec y - \vec X\vec\beta)=\vec y^\top\vec y - 2\vec y^\top\vec X\vec\beta + \vec\beta^\top\vec X^\top\vec X\vec\beta

可以看作是关于 $\vec \beta$ 的二次函数 $F(\vec\beta)$ . 对于确定的常数 $t$ , 函数的等值线 $F(\vec\beta)=t$ 代表了所有使得残差平方和相同的参数 $\vec\beta$ 构成的集合, 这些点关于轨迹便是在 $(\beta_1,\beta_2)$ 平面内的一簇簇椭圆. 越靠近椭圆的中心, 残差平方和越小 . 同时, 对于确定的常数 $s$ , $||\vec\beta||^2=s$ 即可看作是在 $(\beta_1,\beta_2)$ 平面内的圆. 因此如果我们回到带约束条件的优化问题里:

\min_{\vec\beta} || \vec y - \vec X\vec\beta||^2, \quad ||\beta||^2 \leq \lambda^2,

它便代表了在 $(\beta_1,\beta_2)$ 不超出以 $\lambda$ 为半径的圆的前提下使得数据中残差平方和最小. 那么对于确定的 $\lambda$ , 我们很自然地想让 $||\vec\beta||^2=\lambda^2$ 与 $F(\vec\beta)$ 的等值线相切. 切点即为满足条件的最优解. 下面的几幅图很好地阐释了Ridge回归的几何意义:

在实际应用中，正则化参数 $\lambda$ 通常不是由一个简单公式直接给出, 而是通过验证集或交叉验证选择. 在Ridge回归中, 我们一般先有一组候选值 $\lambda_1,\cdots,\lambda_m$ , 常见的参数有 $0.1, 0.05, 0.001$ 等. 我们用这些参数分别拟合 Ridge 模型, 并比较它们在未参与训练的数据上的预测误差. 使预测误差最小的 $\lambda$ 则被视为较合适的正则化参数.

我们还可以从其他角度去解释 $\lambda$ 的不同取值会产生什么样的影响. 回顾第四章第 $4$ 节的知识, 我们知道任何矩阵都存在奇异值分解. 因此对于设计矩阵 $\vec X \in M_{n(m+1)}(\mathbb{R})$ , 存在单位正交矩阵 $\vec U,\vec V$ 与对角矩阵 $\vec\Sigma$ , 使得 $\vec X = \vec U \vec\Sigma \vec V^\top$ , 其中 $\vec U \in M_{n}(\mathbb{R})$ , $\vec V \in M_{(m+1)}(\mathbb{R})$ , $\vec\Sigma \in M_{n(m+1)}(\mathbb{R})$ , 那么

\vec X^\top\vec X =\Big(\vec U\vec\Sigma\vec V^\top\Big)^\top\vec U\vec\Sigma\vec V^\top= \vec V\vec\Sigma^\top\vec U^\top\vec U\vec\Sigma\vec V^\top = \vec V\vec\Sigma^\top\vec\Sigma\vec V^\top.

因此Ridge回归的解 $\hat{\vec\beta}_{R,\lambda} = \Big(\vec X^\top \vec X+\lambda\vec I\Big)^{-1}\vec X^\top\vec y$ 便可以写成

\begin{aligned} \hat{\vec\beta}_{R,\lambda} &= \Big( \vec V\vec\Sigma^\top\vec\Sigma\vec V^\top+\lambda\vec I\Big)^{-1}\cdot\Big(\vec U\vec\Sigma\vec V^\top\Big)^\top\vec y\\ &=\Big[ \vec V\Big( \vec\Sigma^\top\vec\Sigma + \lambda\vec I \Big)\vec V^\top\Big]^{-1}\cdot\vec V\vec\Sigma^\top\vec U^\top\vec y\\ &=\vec V\Big(\vec\Sigma^\top\vec\Sigma+\lambda\vec I\Big)^{-1}\vec\Sigma^\top\vec U^\top\vec y. \end{aligned}

此时响应向量 $\vec y$ 的拟合值即为 $\hat{\vec y} = \vec X \hat{\vec\beta}_{R,\lambda}$ , 即

\hat{\vec y}_{R,\lambda} = \vec U\vec\Sigma\vec V^\top\vec V\Big(\vec\Sigma^\top\vec\Sigma+\lambda\vec I\Big)^{-1}\vec\Sigma^\top\vec U^\top\vec y = \vec U\vec\Sigma\Big(\vec\Sigma^\top\vec\Sigma+\lambda\vec I\Big)^{-1}\vec\Sigma^\top\vec U^\top\vec y.

此时我们设 $\vec X$ 有 $k$ 个非零奇异值 ( $k\leq m+1$ ): $\sigma_1\geq\cdots\geq\sigma_k >0$ . 由于矩阵 $\vec\Sigma,\vec\Sigma^\top$ 为对角矩阵, 那么 $\Big( \vec\Sigma^\top\vec\Sigma + \lambda\vec I \Big)^{-1}$ 也为对角矩阵, 且对角线上的元素为 $\frac{1}{\sigma_i^2+\lambda}$ . 由此可知矩阵 $\vec \Sigma \Big( \vec\Sigma^\top\vec\Sigma + \lambda\vec I \Big)^{-1} \vec\Sigma^\top$ 也为对角矩阵, 其对角线上的元素即为三个对角矩阵对应元素的乘积, 即 $\frac{\sigma_i^2}{\sigma_i^2+\lambda}, i=1,\cdots,k$ . 那么此时我们可以将 $\hat{\vec\beta}$ 写成

\hat{\vec y}_{R,\lambda}= \sum_{i=1}^k \frac{\sigma_i^2}{\sigma_i^2+\lambda}\vec u_i \vec u_i^\top\vec y,\tag{7.4}

其中 $\vec u_i$ 即为矩阵 $\vec U$ 中的列向量, $\frac{\sigma_i^2}{\sigma_i^2+\lambda} \in (0,1)$ . 为了更好地解释(7.4), 我们有必要写出拟合向量在OLS估计量下有着怎样的分解形式. 其实我们不难发现,

\hat{\vec y}_{OLS} = \sum_{i=1}^k \vec u_i \vec u_i^\top \vec y.

因此, Ridge回归可以看作是将OLS中的拟合值进行“缩放”而得到的结果. 在 $\lambda$ 确定时, 缩放的大小则取决于矩阵 $\vec X$ 的奇异值. 奇异值越大, $\frac{\sigma_i^2}{\sigma_i^2+\lambda}$ 一项越接近 $1$ ; 奇异值越小, $\frac{\sigma_i^2}{\sigma_i^2+\lambda}$ 一项越接近 $0$ . 在设计矩阵 $\vec X$ 中, 尤其是当 $(m+1)>>n$ 时, $\vec X$ 的列向量之间会展现出很强的多重共线性 (Multicolinearity), 即可以理解为它们之间彼此线性相关. 在此情况下强烈的多重共线性会使得 $\vec X$ 出现很小的奇异值, 从而导致OLS在这些奇异值方向上的系数极不稳定. Ridge回归的作用可以通俗地理解为: 其更好地保留了 $\vec X$ 中真正对响应变量能起到作用的协变量. 对于一些关联性极小的协变量则削弱它们的权重. 这种从众多协变量中选出部分重要的协变量也被我们称为变量选择 (Variable Selection).

由于Ridge回归具有相对简单的闭合解, 因此我们不妨再研究一下其性质:

定理 7.7

在模型 $\vec y = \vec X\vec\beta+ \vec\epsilon$ 中, 设噪声满足同方差性 $\Var(\epsilon_i) = \sigma^2$ , 设 $\vec X$ 为设计矩阵, $\lambda\geq 0$ , $\hat{\vec\beta}_{Ridge}$ 为参数 $\vec\beta$ 的Rideg回归估计量, 那么

\E[\hat{\vec\beta}_{R,\lambda} | \vec X] = \Big(\vec X^\top\vec X+\lambda\vec I \Big)^{-1}\vec X^\top\vec X\vec\beta, \quad \Var(\hat{\vec\beta}_{R,\lambda}|\vec X) = \sigma^2\Big(\vec X^\top\vec X+\lambda\vec I\Big)^{-1}\vec X^\top\vec X\Big(\vec X^\top\vec X + \lambda\vec I \Big)^{-1}.

该定理读者自证不难. 它告诉我们此时Ridge回归估计量不再为关于参数 $\vec\beta$ 的无偏估计, 但是一个好处便是其方差小于最小二乘估计量的方差.

定理 7.8

对任意的 $\lambda> 0$ , Ridge回归估计量的方差小于最小二乘估计量的方差. 即

\Var(\hat{\vec\beta}_{R,\lambda}|\vec X) \leq \Var(\hat{\vec\beta}_{OLS} | \vec X).

该定理的证明略.

*{LASSO回归}

提起变量选择, 我们就不得不再介绍一种正则化手段: LASSO 回归 (LASSO Regression). LASSO的全称为最小绝对收缩与选择算子(Least Absolute Shrinkage and Selection Operator), 其原理与Ridge回归非常相似, 都是在 $||\vec y-\vec X\vec\beta||^2$ 一项后面加入了一个惩罚项. 根据定义, 通过LASSO回归得到的参数估计 $\hat{\vec\beta}_{LASSO}$ 满足

\hat{\vec\beta}_{LASSO} := \arg\min_{\vec\beta} \left(|| \vec y - \vec X\vec\beta||^2 + \lambda ||\vec\beta||_1\right), \quad\lambda>0,\tag{7.5}

其中我们定义 $||\vec\beta||_1$ 为参数向量 $\vec\beta = \begin{pmatrix} \beta_0 & \beta_1 & \cdots & \beta_m \end{pmatrix}$ 的 $L^1$ 模:

||\vec\beta||_1 := \sum_{i=0}^m |\beta_i|.

对于(7.5)中一般情形的LASSO回归, 我们无法像Ridge回归那样求出一个简单的解. 这是因为在LASSO的惩罚项中有绝对值函数 $|\beta_j|$ , 而该函数在 $\beta_j=0$ 处不可导. 因此我们无法通过简单的求导而得到一个公式. 不过，在一维情形中, 我们可以进行简易分析: 此时(7.5)可以写作

\arg\min_{\beta\in\mathbb R} \left\{ (z-\beta)^2+\lambda|\beta| \right\} = \arg\min_{\beta \in \mathbb{R}} \{ \beta^2 - 2z\beta+\lambda|\beta| \}.\tag{7.6}

方程(7.6)的解为

\hat{\beta}_{LASSO} \begin{cases} z-\frac{\lambda}{2}, & z>\frac{\lambda}{2}\\ 0, & |z|\leq \frac{\lambda}{2}\\ z+\frac{\lambda}{2}, & z<-\frac{\lambda}{2} \end{cases}.

因此，当未正则化估计的 $z$ 的绝对值不超过 $\frac{\lambda}{2}$ 时, LASSO会将对应系数直接压缩为 $0$ . 这样一来便可以减少协变量的数量, 这也正是LASSO可以进行变量选择的原因. 从几何角度来看, 在二维平面 $(\beta_1,\beta_2)$ 中函数 $|\beta_1|+|\beta_2|$ 的等值线即为一个“菱形”. 因此当我们尝试去寻找此时 $|\beta_1|+|\beta_2|$ 与 $||\vec y - \vec X\vec\beta||$ 的等值线相切的点时, 菱形的棱角变更容易成为切点. 此时也便对应了 $\beta_1$ 或 $\beta_2$ 的坐标为零, 也因此达到了变量选择的效果.

{7.3 R Lab}

我们将使用R软件来演示如何运用最小二乘, Ridge与Lasso回归来进行数据分析. 假设我们想要研究一处海滩附近每日的鲨鱼数量. 那么鲨鱼数量便将作为我们的响应变量. 同时假设我们目前考虑的协变量有: 当日的气温, 当日在海滩游泳的人数, 当日海滩附近渔船的数量, 当日AMD股票的均价, 以及当日MLB棒球比赛的场数. 那么不难发现在这些协变量里面最后两个基本上与响应变量没有什么关系, 因此在建立回归模型时我们也便自然希望看到这两个协变量对应的参数几乎为零. 在此, 我们在R软件中进行以下的模拟:

	set.seed(123)
	# 随机生成300个样本
	num <- 300
	#生成协变量
	dat <- data.frame(ships = rnorm(n=num, mean=50, sd=10),
	swimmers = round(rnorm(n=num, mean=500, sd=100)),
	temp = rnorm(n=num, mean=90, sd=2),
	stock_price = runif(n=num, min = 100, max=150),
	mlb_games = runif(n=num, min=0, max=10)
	)

	#响应变量与协变量之间的关系
	sharks <- round(rnorm(n=num, mean = 30, sd=10)+ #噪声项
	-2*dat$ships+0.1*dat$swimmers+1*dat$temp+ 0*dat$stock_price+0*dat$mlb_games)

	dat$sharks <- sharks

	plot(dat)

根据我们假设的模型, 当日AMD股票的均价以及MLB棒球比赛的数量对鲨鱼数量没有任何影响; 每多一艘船便会使得鲨鱼数量下降 $2$ ; 每多 $10$ 位游泳者便会使得鲨鱼数量上升 $1$ ; 气温每提高 $1$ 度也会使得鲨鱼数量上升 $1$ . 最后一行的plot指令可以绘制出在观测数据中不同两个协变量之间的散点图:

我们首先来检验利用最小二乘所得到的线性模型 $\vec y = \vec X\vec\beta + \vec\epsilon$ 的表现. 我们在R中运行如下的代码:

	res <- lm(sharks~., data=dat)

	summary(res)

最后一行的summary指令即可对该模型所得到的参数加以概括:

	Call:
	lm(formula = sharks ~ ., data = dat)

	Residuals:
	Min      1Q  Median      3Q     Max
	-26.346  -6.563   0.620   6.643  33.705

	Coefficients:
	Estimate Std. Error t value Pr(>|t|)
	(Intercept) 24.637398  26.203791   0.940    0.348
	ships       -2.127438   0.060019 -35.446  < 2e-16 ***
	swimmers     0.098607   0.005729  17.212  < 2e-16 ***
	temp         1.114489   0.274962   4.053 6.47e-05 ***
	stock_price  0.023424   0.039430   0.594    0.553
	mlb_games   -0.130108   0.191781  -0.678    0.498
	---
	Signif. codes:
	0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

	Residual standard error: 9.757 on 294 degrees of freedom
	Multiple R-squared:  0.8497,	Adjusted R-squared:  0.8472
	F-statistic: 332.5 on 5 and 294 DF,  p-value: < 2.2e-16

在Estimate一栏中, 我们不难发现对于船只, 游泳人数, 当日气温, AMD股票均价, MLB棒球比赛数量这些协变量参数的估计量分别为 $-2.127$ , $0.099$ , $1.114$ , $0.023$ , $-0.130$ . 这些估计量与真实值已经十分接近. 接下来我们将运行以下的代码进行Ridge回归.

    install.packages("glmnet")
    library(glmnet)

	varmtx <- model.matrix(sharks~.-1, data=dat)
	response <- dat$sharks

	# alpha=0 means ridge regression.
	ridge <- glmnet(scale(varmtx), response, alpha=0)

	# Cross validation to find the optimal lambda penalization
	cv.ridge <- cv.glmnet(varmtx, response, alpha=0)

	# Create a function for labeling the plot below
	lbs_fun <- function(fit, offset_x=1, ...) {
		L <- length(fit$lambda)
		x <- log(fit$lambda[L])+ offset_x
		y <- fit$beta[, L]
		labs <- names(y)
		text(x, y, labels=labs, ...)
	}

	plot(ridge, xvar = "lambda", label=T)
	lbs_fun(ridge)
	abline(v=cv.ridge$lambda.min, col = "red", lty=3)
	abline(v=cv.ridge$lambda.1se, col="blue", lty=3)

该代码会生成如下的参数与正则化参数 $-\log\lambda$ 的关系图像. 从图中我们可以看到当正则化参数趋于零时, 对应的参数估计便趋近于其最小二乘估计; 而随着正则化参数的增大, 对应的参数估计在惩罚项的影响下逐渐收缩, 最后收缩至零.

如果我们运用LASSO回归, 我们可以运行以下的代码:

	lasso <- glmnet(scale(varmtx), response, alpha=1)

	cv.lasso <- cv.glmnet(varmtx, response, alpha=1)

	plot(lasso, xvar = "lambda", label=T)
	lbs_fun(lasso, offset_x = -2)
	abline(v=cv.lasso$lambda.min, col = "red", lty=2)
	abline(v=cv.lasso$lambda.1se, col="blue", lty=2)

通过观察下图, 我们也不难发现在LASSO回归中, 对于任意的 $\lambda$ 取值, 总会有一些协变量的估计为零, 因此和Ridge回归相比LASSO回归更能起到变量选择的作用. 通过观察下图我们也能分辨出哪些协变量是重要的, 而哪些协变量是可有可无的.

值得注意的是, Ridge与LASSO回归同样也适用于设计矩阵满列秩时的情况.