第 15 章 正则化、高维统计与稀疏计算
在传统回归课程中，我们通常假设样本量 \(n\) 明显大于变量个数 \(p\)。但在现代经济统计中，变量数量可能非常多：宏观预测中有大量指标，企业风险评估中有财务、交易和行为变量，问卷调查中也可能包含大量题项。当 \(p\) 较大或变量高度相关时，普通最小二乘估计可能方差很大，甚至无法唯一估计。
正则化的基本思想是在拟合误差之外加入惩罚项，用一点偏差换取更低方差和更稳定预测。本章介绍 Ridge回归和 Lasso 回归，并强调它们的计算实现。
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_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)
}注意，正则化前通常要标准化变量。否则，惩罚项会受到变量量纲影响，尺度大的变量和尺度小的变量被惩罚的程度不可比。
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 法。常用算法之一是坐标下降：每次固定其他系数，只更新一个系数。核心操作是软阈值函数
\[ S(z,\lambda) = \operatorname{sign}(z)(|z|-\lambda)_+, \]
其中 \((a)_+=\max(a,0)\)。
soft_threshold <- function(z, lambda) {
sign(z) * pmax(abs(z) - lambda, 0)
}
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 为什么会产生稀疏解。
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)
lambda_grid <- exp(seq(log(0.02), log(1.5), length.out = 25))
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_best = lambda_best,
selected_variables = sum(abs(fit_lasso$beta) > 1e-8))
#> lambda_best selected_variables
#> 0.1732 4.0000
round(fit_lasso$beta[1:8], 3)
#> [1] 1.164 -0.762 0.773 0.000 0.000 0.000 0.000 0.000由于数据是随机生成的，Lasso 不一定每次都完美选择前三个变量。但整体上，正则化会减少无关变量带来的方差，提高预测稳定性。
15.4 Ridge 与 Lasso 的比较
Ridge 和 Lasso 都会压缩系数，但方式不同：
- Ridge 适合处理许多相关变量共同产生预测信息的情形；
- Lasso 适合希望得到稀疏模型、做变量筛选的情形；
- Ridge 通常不会把系数压缩为 0；
- Lasso 在高度相关变量中可能只选择其中一个，选择结果可能不稳定；
- 两者都需要通过交叉验证或其他方法选择 \(\lambda\)。
下面在同一数据上比较 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 包处理稀疏矩阵。虽然本书不展开稀疏矩阵算法，但要记住一个原则：数据结构会影响算法可行性。高维统计不是只多了几个变量，它会同时改变估计、存储和计算。