岭回归和线性回归的区别 岭回归的参数

基于Cross Validation方法的岭回归参数计算

• 一、Ridge regression
• 二、Cross validation

• 2.1 交叉验证方法介绍
• 2.2 基于LOOCV的岭参数 λ \lambda λ求解

• 三、Generalized cross validation

在利用岭回归(ridge regression)求解回归模型(egression model)时，岭参数的确定可以采用交叉验证(cross validation)方法来确定。很多文献提到交叉验证法求岭参数的目标函数的矩阵表达形式(如下)，但未具体描述如何推导这样的表达，本文主要详细论述该过程，以供参考！

我们先不加说明地描述通过cross validation 求解岭回归参数$\lambda$的最终目标函数的形式：
$P(\lambda)=\frac{1}{n}\|(I-A(\lambda))\boldsymbol{y}\oslash (\boldsymbol{1}-\operatorname{diag}(A(\lambda)))\|^{2}. \tag{*}$

一、Ridge regression

考虑标准回归模型:
$\boldsymbol{y}=X\boldsymbol{\beta}+\boldsymbol{\varepsilon}, \tag{1}$
其中，$X\in \mathbb{R}^{n\times p}$为已知的设计矩阵(design matrix)，$\boldsymbol{y}\in \mathbb{R}^{n}$为响应向量(response vector) ， $\boldsymbol{\varepsilon}\in \mathbb{R}^{n}$ 为 随机白噪声，满足$\mathbb{E}[\boldsymbol{\varepsilon}]=\boldsymbol{0}，\mathbb{E}[\boldsymbol{\varepsilon \varepsilon^{T}}]=\sigma^{2}I$，$\boldsymbol{\beta}\in \mathbb{R}^{p}$未知.

最小二乘估计的目标函数为如下形式：
$J(\boldsymbol{\beta})=\min _{\boldsymbol{\beta}} \|\boldsymbol{y}-X\boldsymbol{\beta}\|^{2}. \tag{2}$
上式关于向量$\boldsymbol{\beta}$求导，并令导数等于$0$，可以求解得到：
$\tilde{\boldsymbol{\beta}}=(X^{T}X)^{-1}X^{T}\boldsymbol{y}. \tag{3}$
当$X^{T}X$为半正定矩阵时，上式中逆无法获得。针对这个问题，引入了岭回归的方法.

岭回归在最小二乘估计的基础上加入了一个扰动项(或$L_{2}$惩罚项, $\lambda\|\boldsymbol{\beta}\|^{2}$). 一方面保证可以求矩阵逆，另一方面控制最优解不会过大. 具体目标函数如下形式：
$J(\boldsymbol{\beta})=\min _{\boldsymbol{\beta}} \|\boldsymbol{y}-X\boldsymbol{\beta}\|^{2}+\lambda \|\boldsymbol{\beta}\|^{2}. \tag{4}$
同样对上式关于$\boldsymbol{\beta}$求导，并令导数等于0，可以求解得到：
$\tilde{\boldsymbol{\beta}}=(X^{T}X+\lambda I)^{-1}X^{T}\boldsymbol{y}. \tag{5}$

说明：

• $L_{2}$惩罚项的加入，使得$X^{T}X+\lambda I$满秩，保证可逆；
• 单位矩阵对角线全为$1$， 像一条山岭一样，这也是岭回归的由来；
• 岭回归方法由于惩罚项的加入使得回归系数$\boldsymbol{\beta}$的估计不是无偏的，因此岭回归方法是一种以降低精度为代价解决病态 (ill-condition) 矩阵问题或者非适定(ill-posed)问题的回归方法.

---

二、Cross validation

首先引入两个概念：

• 模型方差：回归系数$\boldsymbol{\beta}$的方差；
• 模型偏差：$\boldsymbol{\beta}$的预测值与真实值的偏差；

观察(5)式中回归系数$\boldsymbol{\beta}$的估计值与岭参数$\lambda$的取值相关. 当$\lambda$增大时，$|X^{T}X+\lambda I|$就越大，$(X^{T}X+\lambda I)^{-1}$就越小，此时模型的方差就越小；而$\lambda$越大，$\boldsymbol{\beta}$的估计值就越偏离真实值，因此模型的偏差就越大.

因此，岭回归方法的关键是找到一个合理的$\lambda$值来平衡模型的方差和偏差。而针对$\lambda$的估计有大量的文献对其进行了研究，具体可以参考[1]中提到的参考文献。

2.1 交叉验证方法介绍

本文主要介绍交叉验证法。[Wiki]中针对交叉验证方法的介绍包括两大类：

• Exhaustive cross-validation

• Leave-p-out cross-validation （LPOCV）
• Leave-one-out cross-validation (LOOCV) 又称留一交叉验证

• Non-exhaustive cross-validation

• k-fold cross-validation 又称$k$-折交叉验证
• Holdout method
• Repeated random sub-sampling validation

本文使用的方法是留一交叉验证法(LOOCV)，其大体思想：将数据集拆分为$k$组数据，从$k$组数据中挑选出$k-1$组用于模型的训练，剩下一组用于模型的测试。这样一共有$k$个(训练集，测试集)对，而每一对训练集和测试集下都会有对应的模型及模型评分(如均方误差)，进而可以得到一个平均得分。对于$\lambda$的取值，选择平均得分最优的$\lambda$的值。

2.2 基于LOOCV的岭参数$\lambda$求解

David M. Allen在文献[2]中提出求解$\lambda$的Prediction Sum of Squares (PRESS)判据。该判据利用就是LOOCV的思想，具体地讲，就是对于$n$个观测值来说，第$k$个观测值可以用余下的$n-1$个观测值来预测。

文献[1]中有这样的论述：If $\lambda$ is a good choice, then the $k$-th component $[X\beta^{(k)}(\lambda)]_{k}$ of $X\beta^{(k)}(\lambda)$ should be a good predictor of $y_{k}$. 这里$X\beta^{(k)}(\lambda)$ 表示利用除了第$k$个观测值外的$n-1$个观测值得到的$\boldsymbol{y}$的预测值。

因此最小化如下函数即可得到理想的$\lambda$：
$P(\lambda)=\frac{1}{n}\sum_{k=1}^{n}\left([X\beta^{(k)}(\lambda)]_{k}-y_{k}\right)^{2}, \tag{6}$
(注：需要注意的是，变量$[X\beta^{(k)}(\lambda)]$随着$k$的变化而变化，换句话说$[X\beta^{(k)}(\lambda)]_{k}$每次取得是不同向量的不同分量.或者说$\sum_{k=1}^{n}\left([X\beta^{(k)}(\lambda)]_{k}-y_{k}\right)^{2}\neq \|X\beta^{(k)}(\lambda)-\boldsymbol{y}\|^{2}$.)
其中，$\beta^{(k)}(\lambda)$表示忽略第$k$次观测值(也即，将$X$和$\boldsymbol{y}$的第$k$行元素全部置为$0$¹)时利用(5)式计算得到的$\boldsymbol{\beta}$的估计值. 可以用如下形式表达[2] ²：
$\beta^{(k)}(\lambda)=(X^{T}X+\lambda I- \boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T})^{-1}(X^{T}\boldsymbol{y}-y_{i}\boldsymbol{x}_{i}), \tag{7}$
其中，$\boldsymbol{x}_{i}$表示矩阵$X^{T}$的第$i$列, 也可以用$X^{T}\boldsymbol{e}_{i}$来表示.

下面我们证明(6)式可以表述为文章开头(*)中的表达形式。
在正式开始证明之前，先补充一个著名的矩阵逆公式(Sherman-Woodbury-Morrison公式) .

(Sherman-Woodbury-Morrison公式) 已知$A\in \mathbb{R}^{n\times n}$非奇异，$B\in \mathbb{R}^{n\times p}, C\in \mathbb{R}^{p\times n}$, 有下式成立：
$(A+BC)^{-1}=A^{-1}-A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1}. \tag{8}$

令$M:= (X^{T}X+\lambda I)$ 以及 $A(\lambda)=XM^{-1}X^{T}$，考虑如下向量：
$\begin{aligned} X\beta^{(i)}(\lambda)-y&=X(M-\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T})^{-1}(X^{T}\boldsymbol{y}-y_{i}\boldsymbol{x}_{i})-\boldsymbol{y}\\ & = X[M^{-1}+M^{-1}\boldsymbol{x}_{i}(1-\boldsymbol{x}_{i}^{T}M^{-1}\boldsymbol{x}_{i})^{-1}\boldsymbol{x}_{i}^{T}M^{-1}](X^{T}\boldsymbol{y}-y_{i}\boldsymbol{x}_{i})-\boldsymbol{y}\\ &=X\left[M^{-1}+\frac{M^{-1}\boldsymbol{x}_{i}\boldsymbol{x}_{i}^{T}M^{-1}}{(1-\boldsymbol{x}_{i}^{T}M^{-1}\boldsymbol{x}_{i})}\right](X^{T}\boldsymbol{y}-X^{T}\boldsymbol{e}_{i}y_{i})-\boldsymbol{y}\\ &=\left[A(\lambda)+\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}A(\lambda)}{1-[A(\lambda)]_{ii}}\right](\boldsymbol{y}-\boldsymbol{e}_{i}y_{i})-\boldsymbol{y}\\ & = \left[A(\lambda)+\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}A(\lambda)}{1-[A(\lambda)]_{ii}}\right]\boldsymbol{y}-\left[\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}\right]\boldsymbol{y}-\boldsymbol{y}\\ &=\left[A(\lambda)-I+\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}(A(\lambda)-I)\right]\boldsymbol{y}\\ &=\left[\left(I+\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}\right)(A(\lambda)-I)\right]\boldsymbol{y} \end{aligned}$
考虑上述向量的第$i$个分量，也即在上述向量左侧乘$\boldsymbol{e}_{i}^{T}$, 可得：
$\begin{aligned} [X\beta^{(i)}(\lambda)]_{i}-y_{i}=&\left[\boldsymbol{e}_{i}^{T}\left(I+\frac{A(\lambda)\boldsymbol{e}_{i}\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}\right)(A(\lambda)-I)\right]\boldsymbol{y}\\ =&\left[\frac{\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}(A(\lambda)-I)\right]\boldsymbol{y}\\ =&\frac{\boldsymbol{e}_{i}^{T}}{1-[A(\lambda)]_{ii}}(A(\lambda)-I)\boldsymbol{y} \end{aligned}$
因此，
$\left( \begin{array}{c} [X\beta^{(1)}(\lambda)]_{1}-y_{1}\\ \vdots\\ [X\beta^{(n)}(\lambda)]_{n}-y_{n} \end{array} \right)=\left[ \begin{array}{ccc} \frac{1}{1-[A(\lambda)]_{11}} & &\\ & \ddots &\\ & & \frac{1}{1-[A(\lambda)]_{nn}} \end{array} \right](A(\lambda)-I)\boldsymbol{y}$
所以结合(6)式有如下表达式 （类似于文献[3]中3.1.1节中$G(\lambda)$的表达形式）：
$P(\lambda)=\frac{1}{n}\|(I-A(\lambda))\boldsymbol{y}\oslash (\boldsymbol{1}-\operatorname{diag}(A(\lambda)))\|^{2}. \tag{9}$
其中, $\oslash$表示分量除(component-wise division).
或者，等价地写为如下表达形式(类似于文献[1]中公式(2.2)的结果）：
$P(\lambda)=\frac{1}{n}\|B(\lambda)(I-A(\lambda))\boldsymbol{y}\|^{2}, \tag{10}$
其中，
$B(\lambda)=\left[ \begin{array}{ccc} \frac{1}{1-[A(\lambda)]_{11}} & &\\ & \ddots &\\ & & \frac{1}{1-[A(\lambda)]_{nn}} \end{array} \right].$

三、Generalized cross validation

Cross validation在大多数情况是适用的. 但是，一些极端情况是不适用的，如当矩阵$X$为对角矩阵时，此时(6)式中$[X\beta^{(k)}(\lambda)]_{k}$为$0$, 也即此时(6)式变为：
$P(\lambda)=\frac{1}{n}\sum_{k=1}^{n}y_{k}^{2}.$
此时，$P(\lambda)$不存在唯一最小值. 因此，我们可以简单总结为：Cross validation 对于$X$近似对角阵时不适用.

Golub等人在[1]中提出了Generalized cross validation方法. 以下只简单介绍相关思想，具体可参考文献[1]中的介绍.
假设矩阵$X$的SVD分解为:
$X=UDV^{T}.$
假设$W$为可以将循环矩阵变化为对角矩阵的酉矩阵，其中矩阵$W$的$(j,k)$项为：
$[W]_{jk}=\frac{1}{\sqrt{n}}e^{2\pi i j k/n}, \ j,k=1,2,\cdots, n.$
对模型做如下变化：
$\begin{aligned} \tilde{\boldsymbol{y}}=WU^{T}\boldsymbol{y}&=WDV^{T}\boldsymbol{\beta}+WU^{T}\boldsymbol{\varepsilon}\\ &=\tilde{X}\boldsymbol{\beta}+WU^{T}\boldsymbol{\varepsilon}, \end{aligned}$
其中$\tilde{X}=WDV^{T}$.
将$\tilde{y}$和$\tilde{X}$代入(10)式，并记$\tilde{A}(\lambda)$为如下形式：
$\tilde{A}(\lambda)=\tilde{X}(\tilde{X}^{*}\tilde{X}+\lambda I)^{-1}\tilde{X}^{*},$
其中$\tilde{X}^{*}$表示$\tilde{X}$的共轭转置. $\tilde{A}(\lambda)$为循环矩阵且对角元常数，并且$A(\lambda)$和$\tilde{A}(\lambda)$具有相同的特征值.
则$P(\lambda)$变为$V(\lambda)$，具体表达如下：
$V(\lambda)=n\frac{\|(I-\tilde{A}(\lambda))\tilde{\boldsymbol{y}}\|^{2}}{\left(\operatorname{Trace}(I-\tilde{A}(\lambda))\right)^{2}}. \tag{11}$

---

参考文献：
[1] Gene H. Golub, Michael Heath, Grace Wahba. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21:2, 215-223, 1979.
[2] David M. Allen. The relationship between variable selection and data agumentation and a method for prediction. Technomerics, 16, 125-127, 1974.
[3] Bryan J. Kaperick. Diagonal Estimation with probing methods. 2019 [D]

---

1. 文献[1]中错误描述为只忽略$y_{k}$($\ldots$with the $k$th data point $y_{k}$, omitted). ↩︎
2. 文献[2]在第二节Choosing a predictor中描述PRESS时，$\hat{Y}_{i}$的表达式左侧少乘了矩阵$X$. ↩︎