第 15 章 正则化、高维统计与稀疏计算
在传统回归课程中，我们通常假设样本量 \(n\) 明显大于变量个数 \(p\)。但在现代经济统计中，变量数量可能非常多：宏观预测中有大量指标，企业风险评估中有财务、交易和行为变量，问卷调查中也可能包含大量题项。当 \(p\) 较大或变量高度相关时，普通最小二乘估计可能方差很大，甚至无法唯一估计。
正则化的基本思想是在拟合误差之外加入惩罚项，用一点偏差换取更低方差和更稳定预测。本章介绍Ridge 回归和 Lasso 回归，并强调它们的计算实现，尤其是一范数惩罚带来的软阈值函数和 ADMM 思路。
本章有两个层面的目标。第一个目标是统计上的：理解为什么在变量多、相关性强或预测目标优先的场景中，我们愿意牺牲一点训练集拟合程度来换取更稳定的估计。第二个目标是计算上的：理解惩罚项如何改变优化问题的结构。Ridge 的惩罚项光滑，最后仍然落到线性方程组；Lasso 的惩罚项不可导，需要用软阈值、坐标下降或 ADMM 这类更适合不可导优化的工具。
为了让符号简洁，本章默认线性模型为
\[ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, \]
其中 \(\mathbf X\) 是 \(n\times p\) 的设计矩阵。实际建模时通常还要包含截距项。正则化一般不惩罚截距，因为截距只决定整体水平，不代表某个解释变量的影响。一个常用做法是先对 \(\mathbf y\) 中心化，对\(\mathbf X\) 的每一列中心化并标准化，然后只对斜率系数做惩罚；最后再把系数变回原始尺度，并恢复截距。本章代码都采用这个做法。
标准化不仅是技术细节，而是正则化建模的基本要求。若一个变量以“元”为单位，另一个变量以“万元”为单位，同一个 \(\lambda\) 对两个系数的惩罚强度并不公平。标准化后，各变量先被放到可比较的尺度上，惩罚项才主要反映变量与响应变量之间的关系，而不是反映单位选择。
15.1 Ridge 回归
Ridge 回归在最小二乘目标函数中加入二范数惩罚：
\[ \hat{\boldsymbol\beta}_{\text{ridge}} = \arg\min_{\boldsymbol\beta} \left\{ \sum_{i=1}^n(y_i-\mathbf x_i^\top\boldsymbol\beta)^2 +\lambda\sum_{j=1}^p\beta_j^2 \right\}. \]
其中 \(\lambda\geq 0\) 是调节参数。\(\lambda=0\) 时就是普通最小二乘；\(\lambda\) 越大，系数被压缩得越接近 0。若自变量已经中心化，Ridge 解为
\[ \hat{\boldsymbol\beta}_{\text{ridge}} = (\mathbf X^\top\mathbf X+\lambda\mathbf I)^{-1} \mathbf X^\top\mathbf y. \]
这里 \(\mathbf I\) 是单位矩阵。与普通最小二乘相比，\(\lambda\mathbf I\) 改善了矩阵条件，从计算上也提高了稳定性。
从计算角度看，Ridge 有两个重要特点。第一，它保留了二次目标函数的光滑性，因此仍然可以通过解线性方程组得到闭式解。第二，加入 \(\lambda\mathbf I\) 后，即使 \(\mathbf X^\top\mathbf X\) 奇异或接近奇异，\(\mathbf X^\top\mathbf X+\lambda\mathbf I\) 也更容易稳定求解。正因为如此，Ridge 常被看作同时具有统计稳定化和数值稳定化作用的估计方法。
在下面的教学函数中，我们先标准化自变量并中心化响应变量，在标准化尺度上求解 Ridge，再把系数还原到原始尺度。函数返回的 intercept 和 beta 可以直接用于原始数据上的预测。
ridge_fit <- function(X, y, lambda) {
Xs <- scale(X)
yc <- y - mean(y)
center <- attr(Xs, "scaled:center")
scale_x <- attr(Xs, "scaled:scale")
beta_scaled <- solve(crossprod(Xs) + lambda * diag(ncol(Xs)),
crossprod(Xs, yc))
beta <- as.vector(beta_scaled) / scale_x
intercept <- mean(y) - sum(center * beta)
list(intercept = intercept, beta = beta)
}注意，正则化前通常要标准化变量。否则，惩罚项会受到变量量纲影响，尺度大的变量和尺度小的变量被惩罚的程度不可比。
Ridge 的系数通常不会等于 0。即使某个变量的作用很弱，它的系数也只是被连续地压小。这一点适合许多变量共同贡献预测信息的场景，例如一组高度相关的宏观指标共同解释经济增长。但如果研究目标希望得到一个更稀疏、更容易报告的变量集合，Ridge 本身并不能自动完成变量筛选。
15.2 Lasso 回归
Lasso 回归使用一范数惩罚：
\[ \hat{\boldsymbol\beta}_{\text{lasso}} = \arg\min_{\boldsymbol\beta} \left\{ \frac{1}{2n}\sum_{i=1}^n(y_i-\mathbf x_i^\top\boldsymbol\beta)^2 +\lambda\sum_{j=1}^p|\beta_j| \right\}. \]
与 Ridge 不同，Lasso 可以把一些系数精确压缩到 0，因此具有变量选择效果。它适合用于高维预测和初步筛选变量，但不应把 Lasso 选出的变量机械解释为“真实因果变量”。Lasso 的目标函数在\(\beta_j=0\) 处不可导，因此不能直接套用普通 Newton 法。理解它的关键，是先理解软阈值函数。
Lasso 最早由 Tibshirani (1996) 系统提出。名称来自 least absolute shrinkage and selection operator。这里的 absolute 指一范数惩罚，selection 指系数可以被压到 0。它的目标函数可以看成两个力量的平衡：平方损失希望尽量贴近数据，一范数惩罚希望系数总绝对值不要太大。当 \(\lambda\) 增大时，惩罚项变得更重要，模型会从复杂走向稀疏；当 \(\lambda=0\) 时，Lasso 退化为普通最小二乘。
软阈值函数
先看一维问题：
\[ \min_x \left\{\frac{1}{2}(x-z)^2+\lambda |x|\right\}. \]
如果 \(x>0\)，一阶条件给出 \(x=z-\lambda\)，这要求 \(z>\lambda\)；如果 \(x<0\)，一阶条件给出\(x=z+\lambda\)，这要求 \(z<-\lambda\)；如果 \(|z|\leq \lambda\)，最优解落在不可导点 \(x=0\)。因此，这个问题的解可以写成
\[ S(z,\lambda) = \operatorname{sign}(z)(|z|-\lambda)_+, \]
其中 \((a)_+=\max(a,0)\)。这个映射称为软阈值函数（soft-thresholding）。它把绝对值小于\(\lambda\) 的数直接压成 0，把绝对值大于 \(\lambda\) 的数向 0 收缩 \(\lambda\)。从第6.6 节的角度看，软阈值函数正是一范数的近端算子：
\[ \operatorname{prox}_{\lambda \lVert\cdot\rVert_1}(\mathbf z) = S(\mathbf z,\lambda), \]
这里的 \(S(\mathbf z,\lambda)\) 表示对向量的每个分量分别做软阈值。Parikh and Boyd (2014) 对近端算子、软阈值和 Lasso 的关系有系统介绍。
软阈值函数与普通“硬阈值”不同。硬阈值会把小于阈值的数设为 0，却保留大于阈值的数不变；软阈值在保留大系数的同时，还会把它们继续向 0 收缩。因此，Lasso 不仅做变量筛选，也做系数收缩。这种收缩带来偏差，但通常能降低方差，尤其在变量多或共线性强时很有用。
也可以用次梯度来理解这个结果。\(|x|\) 在 \(x=0\) 处的次梯度是区间 \([-1,1]\)。最优性条件为
\[ 0\in x-z+\lambda\partial |x|. \]
当 \(x=0\) 时，这个条件变成 \(z\in \lambda[-1,1]\)，也就是 \(|z|\leq \lambda\)。这正解释了为什么一段区间内的输入会被压成精确的 0。
坐标下降求 Lasso
Lasso 的常用算法之一是坐标下降：每次固定其他系数，只更新一个系数。对第 \(j\) 个坐标，更新问题会退化为上面的一维软阈值问题。因此，坐标下降算法的核心操作就是反复调用软阈值函数。
为了看清更新公式，假设 \(\mathbf X\) 已经标准化、\(\mathbf y\) 已经中心化。固定除 \(\beta_j\) 外的其他系数，记
\[ \mathbf r_j=\mathbf y-\sum_{\ell\neq j}\mathbf x_\ell\beta_\ell \]
为“把第 \(j\) 个变量拿出来以后”的部分残差。此时关于 \(\beta_j\) 的目标函数为
\[ \frac{1}{2n}\lVert \mathbf r_j-\mathbf x_j\beta_j\rVert_2^2 +\lambda|\beta_j|+\text{与 }\beta_j\text{ 无关的项}. \]
展开后可得更新
\[ \beta_j \leftarrow \frac{S(z_j,\lambda)}{\frac{1}{n}\mathbf x_j^\top\mathbf x_j}, \qquad z_j=\frac{1}{n}\mathbf x_j^\top\mathbf r_j. \]
如果每一列标准化到样本方差为 1，则分母接近 1；代码中仍显式保留分母，是为了避免依赖某一种标准化约定。坐标下降算法不断扫描所有坐标，直到系数变化很小。
lasso_cd <- function(X, y, lambda, max_iter = 1000, tol = 1e-6) {
Xs <- scale(X)
yc <- y - mean(y)
center <- attr(Xs, "scaled:center")
scale_x <- attr(Xs, "scaled:scale")
p <- ncol(Xs)
beta <- rep(0, p)
x2 <- colMeans(Xs^2)
for (iter in seq_len(max_iter)) {
beta_old <- beta
for (j in seq_len(p)) {
r_j <- yc - drop(Xs %*% beta) + Xs[, j] * beta[j]
z_j <- mean(Xs[, j] * r_j)
beta[j] <- soft_threshold(z_j, lambda) / x2[j]
}
if (max(abs(beta - beta_old)) < tol) break
}
beta_raw <- beta / scale_x
intercept <- mean(y) - sum(center * beta_raw)
list(intercept = intercept, beta = beta_raw, iter = iter)
}这段代码用于教学，强调算法逻辑。实际项目中应优先使用经过充分测试的专业包，例如 glmnet。但理解坐标下降有助于看清 Lasso 为什么会产生稀疏解。
这段实现为了可读性，每次更新坐标时都重新计算残差，因此并不是最高效版本。成熟软件通常会维护并更新当前残差，避免重复矩阵乘法；还会沿着一串 \(\lambda\) 值从大到小求解，并把上一个 \(\lambda\) 的解作为下一个问题的初始值，这称为 warm start。对 Lasso 来说，warm start 很自然：相邻两个\(\lambda\) 对应的解通常也相近。
ADMM 求 Lasso
ADMM 的思路是把“拟合数据”和“产生稀疏性”拆给两个变量处理。对标准化后的设计矩阵 \(\mathbf X\) 和中心化后的响应变量 \(\mathbf y\)，Lasso 可以改写为
\[ \begin{aligned} \min_{\boldsymbol\beta,\mathbf z}\quad &\frac{1}{2n}\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2 +\lambda\lVert \mathbf z\rVert_1,\\ \text{subject to}\quad &\boldsymbol\beta=\mathbf z. \end{aligned} \]
这里 \(\boldsymbol\beta\) 负责最小二乘拟合，\(\mathbf z\) 负责一范数惩罚，约束\(\boldsymbol\beta=\mathbf z\) 保证最后二者一致。这个改写并没有改变原问题，只是把一个变量复制成两个角色：一个角色适合做二次优化，一个角色适合做软阈值。
ADMM 可以从增广拉格朗日函数出发推导。使用缩放对偶变量 \(\mathbf u\) 后，每一步交替最小化
\[ \frac{1}{2n}\lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2 +\lambda\lVert\mathbf z\rVert_1 +\frac{\rho}{2} \lVert\boldsymbol\beta-\mathbf z+\mathbf u\rVert_2^2, \]
再更新 \(\mathbf u\)。这里的二次项惩罚 \(\boldsymbol\beta\) 和 \(\mathbf z\) 的不一致；\(\mathbf u\) 记录这种不一致的累计方向。由此可得三步更新：
\[ \begin{aligned} \boldsymbol\beta^{k+1} &= \left(\frac{\mathbf X^\top\mathbf X}{n}+\rho\mathbf I\right)^{-1} \left(\frac{\mathbf X^\top\mathbf y}{n}+\rho(\mathbf z^k-\mathbf u^k)\right),\\ \mathbf z^{k+1} &= S\left(\boldsymbol\beta^{k+1}+\mathbf u^k,\frac{\lambda}{\rho}\right),\\ \mathbf u^{k+1} &= \mathbf u^k+\boldsymbol\beta^{k+1}-\mathbf z^{k+1}. \end{aligned} \]
第一行是一个 Ridge 形式的线性方程组；第二行是软阈值；第三行累计两个变量之间的不一致。参数 \(\rho>0\) 不是 Lasso 的调节参数，而是 ADMM 的算法参数。它会影响收敛速度，但在凸问题中通常不改变最终解。
ADMM 的停止条件通常同时看两个残差。原始残差\(\lVert\boldsymbol\beta^{k+1}-\mathbf z^{k+1}\rVert_2\) 衡量约束是否满足；对偶残差\(\rho\lVert\mathbf z^{k+1}-\mathbf z^k\rVert_2\) 衡量软阈值变量是否还在明显变化。教学代码中为了简洁，用两个残差的最大值是否小于 tol 作为停止准则。实际软件会根据变量维度、绝对误差和相对误差设置更细的阈值。
lasso_admm <- function(X, y, lambda, rho = 1,
max_iter = 2000, tol = 1e-5) {
Xs <- scale(X)
yc <- y - mean(y)
center <- attr(Xs, "scaled:center")
scale_x <- attr(Xs, "scaled:scale")
n <- nrow(Xs)
p <- ncol(Xs)
beta <- rep(0, p)
z <- rep(0, p)
u <- rep(0, p)
A <- crossprod(Xs) / n + rho * diag(p)
b <- as.vector(crossprod(Xs, yc) / n)
chol_A <- chol(A)
solve_A <- function(rhs) {
backsolve(chol_A, forwardsolve(t(chol_A), rhs))
}
for (iter in seq_len(max_iter)) {
z_old <- z
beta <- as.vector(solve_A(b + rho * (z - u)))
z <- soft_threshold(beta + u, lambda / rho)
u <- u + beta - z
primal_resid <- sqrt(sum((beta - z)^2))
dual_resid <- rho * sqrt(sum((z - z_old)^2))
if (max(primal_resid, dual_resid) < tol) break
}
beta_raw <- z / scale_x
intercept <- mean(y) - sum(center * beta_raw)
list(intercept = intercept, beta = beta_raw, iter = iter,
primal_resid = primal_resid, dual_resid = dual_resid)
}对小型教学例子，坐标下降代码更直观；对更复杂的模型，ADMM 的优势在于可以把不同结构拆开处理，并且有时可以缓存线性方程组的分解或并行计算各个子问题。
坐标下降和 ADMM 都可以求解 Lasso，但它们强调的问题结构不同。坐标下降把注意力放在“逐个系数更新”上，每一步都把一个高维问题变成一维软阈值问题，因此特别适合说明不可导惩罚为什么会产生稀疏解。ADMM 则把注意力放在“变量分裂”上，每一轮在线性方程组和软阈值之间交替，并用残差检查两个变量是否已经足够一致，因此特别适合说明约束、近端算子和分块优化之间的关系。
如果只是普通 Lasso，坐标下降通常非常有效；如果模型还包含其他约束、分组结构、矩阵变量或需要分布式计算，ADMM 的变量分裂思想就更有价值。
15.3 调节参数的选择
正则化强度 \(\lambda\) 控制偏差与方差之间的权衡。常用方法是交叉验证：在一组候选 \(\lambda\) 上分别拟合模型，选择验证误差最小的值。
下面构造一个模拟宏观预测数据集。变量数较多，但只有少数变量真正影响目标变量。为了模拟变量相关性，所有变量都叠加了同一个公共因子。
set.seed(2026)
n <- 140
p <- 30
factor_common <- rnorm(n)
X <- matrix(rnorm(n * p), n, p) + 0.4 * factor_common
beta_true <- c(1.2, -1.0, 0.8, rep(0, p - 3))
y <- drop(X %*% beta_true) + rnorm(n, sd = 1.5)用 5 折交叉验证选择 Lasso 的 \(\lambda\)：
set.seed(1)
Xs_full <- scale(X)
yc_full <- y - mean(y)
lambda_max <- max(abs(crossprod(Xs_full, yc_full) / n))
lambda_grid <- exp(seq(log(lambda_max), log(lambda_max * 0.02),
length.out = 30))
K <- 5
fold <- sample(rep(seq_len(K), length.out = n))
cv_mse <- numeric(length(lambda_grid))
for (l in seq_along(lambda_grid)) {
loss <- numeric(K)
for (k in seq_len(K)) {
train <- fold != k
valid <- fold == k
fit <- lasso_cd(X[train, ], y[train], lambda = lambda_grid[l])
pred <- fit$intercept + X[valid, ] %*% fit$beta
loss[k] <- mean((y[valid] - pred)^2)
}
cv_mse[l] <- mean(loss)
}
lambda_best <- lambda_grid[which.min(cv_mse)]
fit_lasso <- lasso_cd(X, y, lambda = lambda_best)
c(lambda_max = lambda_max,
lambda_best = lambda_best,
selected_variables = sum(abs(fit_lasso$beta) > 1e-8))
#> lambda_max lambda_best selected_variables
#> 1.1525 0.1744 4.0000
round(fit_lasso$beta[1:8], 3)
#> [1] 1.162 -0.760 0.772 0.000 0.000 0.000 0.000 0.000这里的 lambda_max 是一个有用的参考值：当 \(\lambda\geq\lambda_{\max}\) 时，Lasso 的零向量满足最优性条件，因此所有斜率系数都会被压到 0。实际搜索 \(\lambda\) 时，常从 lambda_max 开始逐渐减小，这样能覆盖从“完全稀疏”到“较复杂模型”的路径。
由于数据是随机生成的，Lasso 不一定每次都完美选择前三个变量。但整体上，正则化会减少无关变量带来的方差，提高预测稳定性。这里输出前 8 个系数，是为了观察前三个真实非零变量是否被保留，以及其他变量是否被压到 0。
用同一个 \(\lambda\) 运行 ADMM，可以看到它给出的系数与坐标下降基本一致：
fit_admm <- lasso_admm(X, y, lambda = lambda_best, rho = 1)
c(cd_iter = fit_lasso$iter,
admm_iter = fit_admm$iter,
max_diff = max(abs(fit_lasso$beta - fit_admm$beta)),
admm_nonzero = sum(abs(fit_admm$beta) > 1e-8))
#> cd_iter admm_iter max_diff admm_nonzero
#> 8.0e+00 3.1e+01 7.1e-06 4.0e+00坐标下降和 ADMM 是两个不同的算法，迭代次数不能直接比较，因为“一次迭代”的含义不同。这里更重要的检查是 max_diff：如果两个算法实现正确，并且停止阈值足够小，它们得到的系数应该非常接近。
15.4 Ridge 与 Lasso 的比较
Ridge 和 Lasso 都会压缩系数，但方式不同：
- Ridge 适合处理许多相关变量共同产生预测信息的情形；
- Lasso 适合希望得到稀疏模型、做变量筛选的情形；
- Ridge 通常不会把系数压缩为 0；
- Lasso 在高度相关变量中可能只选择其中一个，选择结果可能不稳定；
- 两者都需要通过交叉验证或其他方法选择 \(\lambda\)。
在解释结果时，还需要注意以下几点。第一，正则化后的系数有偏，不能像普通最小二乘那样直接套用常规标准误和 \(t\) 检验。第二，Lasso 的变量选择受样本扰动影响，尤其在多个变量高度相关时，算法可能在一组相似变量中选择其中一个。第三，预测任务中更关心测试误差；解释任务中则需要结合研究设计、变量含义和稳健性分析。
下面在同一数据上比较 Ridge 和 Lasso 的训练误差。这里只是演示，不代表泛化误差。
fit_ridge <- ridge_fit(X, y, lambda = 20)
pred_ridge <- fit_ridge$intercept + X %*% fit_ridge$beta
pred_lasso <- fit_lasso$intercept + X %*% fit_lasso$beta
pred_ols <- fitted(lm(y ~ X))
c(OLS = mean((y - pred_ols)^2),
Ridge = mean((y - pred_ridge)^2),
Lasso = mean((y - pred_lasso)^2))
#> OLS Ridge Lasso
#> 1.932 2.029 2.331训练误差最小的往往是最灵活的模型，但这不等于预测最好。因此，本章应与第 14章的交叉验证一起理解。
15.5 稀疏计算简介
高维数据不仅带来统计问题，也带来计算问题。若一个矩阵大部分元素为 0，称为稀疏矩阵。许多经济和管理数据都可能具有稀疏性，例如用户–商品购买矩阵、地区–行业投入产出矩阵、企业–银行借贷网络。
稀疏矩阵不应按普通密集矩阵存储，否则会浪费大量内存。R 中常用 Matrix 包处理稀疏矩阵。虽然本书不展开稀疏矩阵算法，但要记住一个原则：数据结构会影响算法可行性。高维统计不是只多了几个变量，它会同时改变估计、存储和计算。