Ridge Regression

内容提要

一周以前（3月7日）考计量期中，第一题是证明岭回归估计量的渐进分布的，我没有完整地写出来，这里面还是有些技巧，之前我并不熟悉，现在总结在这里。

岭回归简介

Ridge regression，中文译名「岭回归」，是一种利用 正则化方法 处理 不适定问题 的技术。我们知道，使用 OLS 时是要求自变量之间不存在完全多重共线性的，即回归元矩阵必须是列满秩的，否则最小二乘问题没有唯一解；但如果某些变量间的相关性很高，甚至接近完全相关，OLS 估计量的方差将会相当大，估计精度大大下降。这个问题在待估参数很多时常常发生。岭回归通过对系数向量的 $L^2$ 范数（的平方）施加一个惩罚，控制系数的整体大小，达到以一定程度的偏误换取更低方差的目的（见 bias-variance tradeoff）。

Formulation

考虑线性模型：$y_i = \bm{x}_i^\prime\bm{\beta} + \varepsilon_i$，$i=1,\dots,N$。其中 $\bm{x}_i$ 和 $\bm{\beta}$ 都是 $p$ 维向量。我们更喜欢表示为矩阵形式：

$$ \bm{Y} = \bm{X}\bm{\beta} + \bm{\mathcal{E}} $$

其中 $\bm{Y} = (y_1,\dots,y_N)^\prime$，$\bm{X} = (\bm{x}_1,\dots,\bm{x}_N)^\prime$，$\bm{\mathcal{E}} = (\varepsilon_1,\dots,\varepsilon_N)^\prime$。假设 $\mathbb{E}(\varepsilon_i|\bm{x}_i) = 0$ 以及 $\mathbb{E}(\varepsilon_i^2|\bm{x}_i) = \sigma^2$，并假设 $\{(\bm{x}_i,y_i)\}_{i=1}^n$ 是独立同分布的。

岭回归估计量 $\hat{\bm{\beta}}$ 来自如下带有惩罚项的最小二乘问题：

$$ \min_{\bm{\beta}}\quad \frac1n(\bm{Y}-\bm{X}\bm{\beta})^\prime(\bm{Y}-\bm{X}\bm{\beta}) + \lambda\bm{\beta}^\prime\bm{\beta} \tag{$*$} $$

这里 $\lambda$ 是一个参数，它衡量了惩罚的大小。当 $\lambda=0$ 时，该问题就是 OLS；当 $\lambda\to\infty$ 时，$\bm{\beta} = \bm{0}$。事实上，这个问题等价于一个有约束的最小二乘：

$$ \min_{\bm{\beta}}\quad (\bm{Y}-\bm{X}\bm{\beta})^\prime(\bm{Y}-\bm{X}\bm{\beta}) \qquad \text{s.t.}\quad \bm{\beta}^\prime\bm{\beta} = c^2 $$

其中 $c$ 为某个非负常数。可见我们是在 $p$ 维 Euclidean 空间中半径为 $c$ 的球面上寻找最小二乘估计量。

由 $(*)$，一阶条件为 $n^{-1}(-2\bm{X}^\prime\bm{Y}+2\bm{X}^\prime\bm{X}\bm{\beta})+2\lambda\bm{\beta} = \bm{0}$，于是岭回归估计量

$$ \hat{\bm{\beta}} = \left(\frac1n\bm{X}^\prime\bm{X}+\lambda\bm{I}_p\right)^{-1}\frac1n\bm{X}^\prime\bm{Y} $$

其中 $\bm{I}_p$ 是 $p$ 阶单位阵。

渐进分布

我们把 $\bm{Y} = \bm{X}\bm{\beta} + \bm{\mathcal{E}}$ 带入 $\hat{\bm{\beta}}$ 的表达式，并做适当整理变形，得到

$$ \hat{\bm{\beta}} - \bm{\beta} = \left(\frac1n\bm{X}^\prime\bm{X}+\lambda\bm{I}_p\right)^{-1}\frac1n\bm{X}^\prime\bm{\mathcal{E}} - \left(\frac1n\bm{X}^\prime\bm{X}+\lambda\bm{I}_p\right)^{-1}\lambda\bm{\beta} $$

为简化符号，记 $\bm{B}_i = \bm{x}_i\bm{x}_i^\prime+\lambda\bm{I}_p$，$\bar{\bm{B}}_n = n^{-1}\sum_{i=1}^n\bm{B}_i$，其无条件期望 $\mathbb{E}(\bm{x}_i\bm{x}_i^\prime)+\lambda\bm{I}_p$ 记作 $\bm{B}$；此外，记 $\bm{Q} = \mathbb{E}(\bm{x}_i\bm{x}_i^\prime)$。

为应用中心极限定理（CLT），需要将上式整理成某个样本均值的形式。事实上，我们应当考虑构造一个更大的随机向量，导出其渐进分布，然后应用 Cramér–Wold device 得到估计量的渐进分布。如果我们构造 $(\bm{x}_i^\prime\varepsilon_i,\lambda\bm{\beta}^\prime)^\prime$，则无法直接应用 CLT，因为该向量均值非零（除非 $\bm{\beta} = \bm{0}$）。正确的方式是作如下变形：

$$ \begin{aligned} \hat{\bm{\beta}} - \bm{\beta} + \lambda\bm{B}^{-1}\bm{\beta} &= \bar{\bm{B}}_n^{-1}\left[\frac1n\bm{X}^\prime\bm{\mathcal{E}}-(\bm{I}_p-\bar{\bm{B}}_n\bm{B}^{-1})\lambda\bm{\beta}\right] \\ &= \bar{\bm{B}}_n^{-1}\frac1n\sum_{i=1}^n\Bigl[\bm{x}_i\varepsilon_i-(\bm{I}_p-\bm{B}_i\bm{B}^{-1})\lambda\bm{\beta}\Bigr] \end{aligned} $$

Note

我们可以将 summand 视为一体，计算其均值和方差，然后直接应用 CLT，这本质上和 Cramér–Wold device 相同。

考虑随机向量 $(\bm{c}_i^\prime,\bm{e}_i^\prime)^\prime$，其中 $\bm{c}_i = \bm{x}_i\varepsilon_i$，$\bm{e}_i = (\bm{I}_p-\bm{B}_i\bm{B}^{-1})\lambda\bm{\beta}$。易见 $\{(\bm{c}_i^\prime,\bm{e}_i^\prime)^\prime\}_{i=1}^n$ 是独立同分布的，均值为 $\bm{0}$，协方差矩阵

$$ \bm{V} = \begin{pmatrix} \sigma^2\bm{Q} & \bm{0} \\ \bm{0} & \lambda^2\bm{D} \end{pmatrix} $$

其中 $\bm{D} = \mathbb{E}(\bm{B}_i\bm{B}^{-1}\bm{\beta\beta}^\prime\bm{B}^{-1}\bm{B}_i)-\bm{\beta}\bm{\beta}^\prime$。由 CLT，有

$$ \frac1{\sqrt{n}}\sum_{i=1}^n\begin{pmatrix} \bm{c}_i \\ \bm{e}_i \end{pmatrix} \to_d N(\bm{0},\bm{V}) $$

于是，利用 Cramér–Wold device，有

$$ \begin{aligned} \frac1{\sqrt{n}}\sum_{i=1}^n\Bigl[\bm{x}_i\varepsilon_i-(\bm{I}_p-\bm{B}_i\bm{B}^{-1})\lambda\bm{\beta}\Bigr] &= (1,-1)\frac1{\sqrt{n}}\sum_{i=1}^n\begin{pmatrix} \bm{c}_i \\ \bm{e}_i \end{pmatrix} \begin{pmatrix} 1 \\ -1 \end{pmatrix} \\ &\to_d N(\bm{0}, \sigma^2\bm{Q}+\lambda^2\bm{D}) \end{aligned} $$

最后，$\bar{\bm{B}}_n\to_p\bm{B}$，根据 Slutsky 定理，有

$$ \sqrt{n}(\hat{\bm{\beta}} - \bm{\beta} + \lambda\bm{B}^{-1}\bm{\beta}) \to_d N(\bm{0},\bm{\Omega}) $$

其中 $\bm{\Omega} = \bm{B}^{-1}(\sigma^2\bm{Q}+\lambda^2\bm{D})\bm{B}^{-1}$。

Note

值得一提，Knight and Fu (2000) 也给出了岭回归估计量的渐进分布，但其做法根本上是不同的，且不适用于我们这里的情形。由于他们并未在 $(*)$ 中第一项前添加 $1/n$，所以按照其文中定义，$\lambda_n = n\lambda$ 并非 $O(\sqrt{n})$，故其 Theorem 2 不能适用。

延伸思考

如果我们想构造渐进方差 $\bm{\Omega}$ 的一致估计，则需要构造 $\bm{B},\bm{Q},\bm{D}$ 的一致估计，前两者容易，而 $\bm{D}$ 较为困难。

记 $\bm{A} = \bm{B}^{-1}\bm{\beta\beta}^\prime\bm{B}^{-1}$。
$\bm{\beta}$ 的一致估计 $\bar{\bm{\beta}} = (\bm{I}_p - \lambda\bar{\bm{B}}_n^{-1})^{-1}\hat{\bm{\beta}}$。
$\bm{A}$ 的一致估计 $\bar{\bm{A}} = \bar{\bm{B}}_n^{-1}\bar{\bm{\beta}}\bar{\bm{\beta}}^\prime\bar{\bm{B}}_n^{-1}$。

乍一看，似乎不能用 $n^{-1}\sum_i\bm{B}_i\bar{\bm{A}}\bm{B}_i$ 作为 $\mathbb{E}(\bm{B}_i\bm{A}\bm{B}_i)$ 的一致估计，因为 $\bar{\bm{A}}$ 本身就是一个估计量，$\bm{B}_i\bar{\bm{A}}\bm{B}_i$ 将不是一列独立随机矩阵，也就无法引用大数定律得到一致性。但实际上是可行的，只需利用 trace 的性质对每个 entry 进行分析：

$\bm{B}_i\bar{\bm{A}}\bm{B}_i$ 的第 $(k,\ell)$ 元素是 $\bm{b}_{ik}^\prime\bar{\bm{A}}\bm{b}_{i\ell}$，$\bm{b}_{ik}^{\prime}$ 是 $\bm{B}_i$ 的第 $k$ 行，注意 $\bm{B}_i$ 是对称矩阵。
$\frac1n\sum_i\text{tr}(\bm{b}_{ik}^\prime\bar{\bm{A}}\bm{b}_{i\ell}) = \text{tr}(\frac1n\sum_i\bm{b}_{i\ell}\bm{b}_{ik}^\prime\bar{\bm{A}})$ 依概率收敛到 $\mathbb{E}(\bm{B}_i\bm{A}\bm{B}_i)$ 的第 $(k,\ell)$ 元素，也就是 $\mathbb{E}(\bm{b}_{ik}^\prime\bm{A}\bm{b}_{i\ell}) = \text{tr}(\mathbb{E}[\bm{b}_{i\ell}\bm{b}_{ik}^\prime]\bm{A})$。

这和 OLS 中异方差时对渐进方差的估计有些相似，比如都是三明治形式。

参考

“Ridge regression.” Wikipedia, accessed March 17, 2023.
van Wieringen, W. N. “Lecture Notes on Ridge Regression.” arXiv preprint, submitted in 2015. arXiv:1509.09169.
Knight, Keith and Wenjiang Fu. 2000. “Asymptotics for Lasso-type Estimators.” The Annals of Statistics, 28 (5): 1356-1378.
Hansen, Bruce. 2022. Econometrics. Princeton University Press.

最后修改于 2024-09-04