第 14 章 正则化、高维统计与稀疏计算

当解释变量较多、变量之间高度相关，或者参数个数接近甚至超过样本量时，普通最小二乘会遇到两类困难。统计上，系数可能随样本扰动剧烈变化，样本内拟合改善却不能转化为预测改进；计算上，设计矩阵可能秩亏或接近秩亏，线性方程组难以稳定求解。

正则化在拟合损失之外加入对系数规模或结构的约束。Ridge 连续压缩系数，Lasso 可以把部分系数压到0，Elastic Net 则把两类惩罚结合起来。惩罚带来偏差，同时可能显著降低方差和预测误差。这种权衡必须由样本外评价决定，而不能只看训练误差。

本章关注从统计目标到计算实现的完整链条：统一目标函数的尺度，解释标准化与截距处理，利用 SVD 计算Ridge，利用软阈值和坐标下降计算 Lasso，并用最优性条件检查算法是否真正收敛。随后按照第13 章的原则选择调节参数，分析系数与选择稳定性，并讨论稀疏矩阵怎样改变存储和计算方式。

14.1 学习目标

完成本章后，读者应当能够：

  1. 从偏差—方差权衡解释正则化为什么可能改善预测；
  2. 说明截距不受惩罚以及自变量标准化的原因；
  3. 利用 SVD 推导并稳定计算 Ridge 解和有效自由度；
  4. 从次梯度条件推导软阈值和 Lasso 坐标更新；
  5. 使用 KKT 条件检查 Lasso 或 Elastic Net 的计算结果；
  6. 解释 Elastic Net 在相关变量和高维问题中的作用；
  7. 使用交叉验证选择惩罚强度，并避免预处理与调参泄漏；
  8. 区分统计稀疏性与矩阵存储稀疏性；
  9. 报告参数选择、收敛诊断与选择稳定性的关键证据。

14.2 正则化问题的建立

正则化首先是一种统计取舍，而不只是让矩阵变得可逆的数值技巧。设真实回归函数为 \(f(\mathbf x)\)，训练数据\(\mathcal D\) 产生估计量 \(\hat f_{\mathcal D}\)。在平方损失下，对固定预测点 \(\mathbf x\)，若新观测噪声方差为\(\sigma^2(\mathbf x)\)，则

\[ \operatorname{E}_{\mathcal D,Y^\ast}\!\left[ \{Y^\ast-\hat f_{\mathcal D}(\mathbf x)\}^2\mid \mathbf X^\ast=\mathbf x \right] =\sigma^2(\mathbf x) +\left\{\operatorname{E}_{\mathcal D}[\hat f_{\mathcal D}(\mathbf x)]-f(\mathbf x)\right\}^2 +\operatorname{Var}_{\mathcal D}\{\hat f_{\mathcal D}(\mathbf x)\}. \]

三项依次是不可约噪声、偏差平方和估计方差。弱惩罚通常偏差较小，却可能对样本扰动十分敏感；增强惩罚会引入收缩偏差，同时降低估计方差。预测误差是否下降取决于两者的合计变化，因此调节参数必须通过样本外预测来选择。对于解释性分析，还要另外考察系数含义和选择不确定性。

14.2.1 统一目标函数

设中心化后的响应为 \(\mathbf y_{\mathrm{c}}\in\mathbb R^n\)，标准化后的设计矩阵为\(\mathbf X_{\mathrm{s}}\in\mathbb R^{n\times p}\)，斜率参数向量为\(\boldsymbol\beta\in\mathbb R^p\)。本章统一使用目标函数

\[\begin{equation} Q_{\lambda,\alpha}(\boldsymbol\beta) =\frac{1}{2n} \lVert\mathbf y_{\mathrm{c}}-\mathbf X_{\mathrm{s}}\boldsymbol\beta\rVert_2^2 +\lambda\left\{ \alpha\lVert\boldsymbol\beta\rVert_1 +\frac{1-\alpha}{2}\lVert\boldsymbol\beta\rVert_2^2 \right\}, \tag{14.1} \end{equation}\]

其中 \(\lambda\geq 0\) 控制总惩罚强度，\(0\leq\alpha\leq1\) 控制两种惩罚的比例。

\(\alpha\) 方法 主要性质
\(0\) Ridge 二次惩罚，连续压缩，通常不产生零系数
\(1\) Lasso 一范数惩罚，可产生零系数
\(0<\alpha<1\) Elastic Net 同时具有稀疏化和二次稳定化作用

目标函数中的 \(1/(2n)\) 很重要。有的软件使用 \(\operatorname{RSS}\)，有的使用\(\operatorname{RSS}/2\) 或 \(\operatorname{RSS}/(2n)\)；同一个数值的 \(\lambda\) 在这些约定下代表不同惩罚强度。本章所有推导和代码都采用式 (14.1) 的尺度。

14.2.2 截距、标准化与原尺度还原

原始尺度上的模型写作

\[ Y_i=\beta_0+\mathbf x_i^\top\boldsymbol\beta+\varepsilon_i. \]

其中随机响应 \(Y_i\) 的观测值记为 \(y_i\)。截距描述响应的整体水平。若惩罚截距，简单改变响应变量的计量原点也会改变斜率估计。常用做法是令

\[ y_{\mathrm{c},i}=y_i-\bar y, \qquad x_{\mathrm{s},ij}=\frac{x_{ij}-\bar x_j}{s_j}, \]

只对标准化尺度上的斜率系数进行惩罚。式 (14.1) 中的\(\boldsymbol\beta\) 始终表示标准化尺度的斜率参数；其估计向量记为\(\widehat{\boldsymbol\beta}\)，第 \(j\) 个分量为 \(\hat\beta_j\)。为与之区分，把还原到原始变量尺度的斜率和截距分别记为 \(\hat\beta_{\mathrm{raw},j}\) 与 \(\hat\beta_{\mathrm{raw},0}\)：

\[ \hat\beta_{\mathrm{raw},j}=\frac{\hat\beta_j}{s_j}, \qquad \hat\beta_{\mathrm{raw},0} =\bar y-\sum_{j=1}^p\bar x_j\hat\beta_{\mathrm{raw},j}. \]

标准化使不同量纲的变量受到可比较的惩罚。若一列变量以“元”计量，另一列以“万元”计量，不经标准化时，相同的实质效应会对应不同大小的系数，惩罚结果将随单位选择改变。二元指示变量是否标准化要结合应用解释和软件约定决定；无论采用哪一种方案，都必须在验证时保持一致。

下面的辅助函数使用均方根尺度\(s_j=\{n^{-1}\sum_{i=1}^n(x_{ij}-\bar x_j)^2\}^{1/2}\)，因此标准化后每列满足\(n^{-1}\sum_{i=1}^n x_{\mathrm{s},ij}^2=1\)。R 的 scale() 默认采用分母 \(n-1\)，两种约定都可以，但推导和代码不能混用。

standardize_xy <- function(X, y) {
  X <- as.matrix(X)
  storage.mode(X) <- "double"
  y <- as.numeric(y)

  if (nrow(X) != length(y) || nrow(X) < 2L || ncol(X) < 1L) {
    stop("X、y 的维数不符合要求。")
  }
  if (any(!is.finite(X)) || any(!is.finite(y))) {
    stop("X 和 y 不能包含缺失值或非有限值。")
  }

  x_center <- colMeans(X)
  X_centered <- sweep(X, 2, x_center, FUN = "-")
  x_scale <- sqrt(colMeans(X_centered^2))
  if (any(!is.finite(x_scale)) || any(x_scale <= 0)) {
    stop("X 含有常数列；请在正则化前删除或单独处理。")
  }

  y_center <- mean(y)
  list(
    X = sweep(X_centered, 2, x_scale, FUN = "/"),
    y = y - y_center,
    x_center = x_center,
    x_scale = x_scale,
    y_center = y_center
  )
}

restore_coefficients <- function(beta_standardized, prep) {
  beta <- as.numeric(beta_standardized) / prep$x_scale
  intercept <- prep$y_center - sum(prep$x_center * beta)
  list(intercept = intercept, beta = beta)
}

predict_penalized <- function(object, newX) {
  newX <- as.matrix(newX)
  drop(object$intercept + newX %*% object$beta)
}

在交叉验证中，x_center、x_scale 和 y_center 都必须只用当前训练折估计，再用于验证折。先对全样本标准化再分折，会把验证数据的信息带入拟合过程。

14.2.3 贯穿案例：相关的高维预测变量

下面生成训练数据和独立测试数据。前 30 个变量分成两个相关组，真实信号只来自前 5 个变量，其余变量不影响条件均值。训练样本量为 140，变量数为 60；这种设置既包含共线性，也为无关变量吸收训练噪声留下空间。

simulate_regularized_data <- function(n, p = 60L) {
  if (p < 30L) stop("这个模拟要求 p 至少为 30。")

  X <- matrix(rnorm(n * p), nrow = n, ncol = p)
  common_1 <- rnorm(n)
  common_2 <- rnorm(n)
  X[, 1:15] <- X[, 1:15] + 0.8 * common_1
  X[, 16:30] <- X[, 16:30] + 0.8 * common_2
  colnames(X) <- paste0("x", seq_len(p))

  beta_true <- c(1.2, -1.0, 0.8, 0.6, -0.5, rep(0, p - 5L))
  y <- drop(X %*% beta_true) + rnorm(n, sd = 1.5)

  list(X = X, y = y, beta = beta_true)
}

set.seed(2027)
regularized_train <- simulate_regularized_data(140)
regularized_test <- simulate_regularized_data(3000)

X_train <- regularized_train$X
y_train <- regularized_train$y
X_test <- regularized_test$X
y_test <- regularized_test$y
beta_true <- regularized_train$beta

测试数据在调节参数选择结束之前不会使用。因为这是模拟实验，真实系数可用于教学核对；实际数据分析通常没有这项信息。

14.3 Ridge 回归

14.3.1 解的结构与 SVD

令式 (14.1) 中 \(\alpha=0\)，Ridge 估计为 (Hoerl and Kennard 1970)

\[ \widehat{\boldsymbol\beta}_{\mathrm{ridge}}(\lambda) =\arg\min_{\boldsymbol\beta} \left\{ \frac{1}{2n}\lVert\mathbf y_{\mathrm{c}}-\mathbf X_{\mathrm{s}}\boldsymbol\beta\rVert_2^2 +\frac{\lambda}{2}\lVert\boldsymbol\beta\rVert_2^2 \right\}. \]

一阶条件给出

\[ \left(\frac{\mathbf X_{\mathrm{s}}^\top\mathbf X_{\mathrm{s}}}{n} +\lambda\mathbf I_p\right)\widehat{\boldsymbol\beta}_{\mathrm{ridge}}(\lambda) =\frac{\mathbf X_{\mathrm{s}}^\top\mathbf y_{\mathrm{c}}}{n}. \]

这个方程说明二次惩罚会把特征值整体向上平移，但教学代码不必据此显式构造\(\mathbf X_{\mathrm{s}}^\top\mathbf X_{\mathrm{s}}\)。令\(r=\operatorname{rank}(\mathbf X_{\mathrm{s}})\)，其紧致 SVD 为

\[ \mathbf X_{\mathrm{s}}=\mathbf U\mathbf D\mathbf V^\top, \qquad \mathbf D=\operatorname{diag}(d_1,\ldots,d_r), \qquad d_j>0,\quad j=1,\ldots,r. \]

其中 \(\mathbf U\in\mathbb R^{n\times r}\)、\(\mathbf V\in\mathbb R^{p\times r}\) 的列分别标准正交，即

\[ \mathbf U^\top\mathbf U=\mathbf I_r, \qquad \mathbf V^\top\mathbf V=\mathbf I_r, \]

Ridge 解可写为

\[ \widehat{\boldsymbol\beta}_{\mathrm{ridge}}(\lambda) =\mathbf V \operatorname{diag}\left\{ \frac{d_j}{d_j^2+n\lambda}:j=1,\ldots,r \right\} \mathbf U^\top\mathbf y_{\mathrm{c}}. \]

普通最小二乘在小奇异值方向上含有 \(1/d_j\)，会放大数据扰动；Ridge 把它替换为\(d_j/(d_j^2+n\lambda)\)。当 \(d_j\) 很小时，后者趋近于 0，因此不稳定方向受到更强压缩。即使\(p>n\) 或设计矩阵秩亏，只要 \(\lambda>0\)，Ridge 解仍然唯一。

Ridge 拟合值对 \(\mathbf y_{\mathrm{c}}\) 是线性的，其有效自由度为

\[ \operatorname{df}_{\mathrm{ridge}}(\lambda) =\sum_{j=1}^r\frac{d_j^2}{d_j^2+n\lambda}. \]

\(\lambda=0\) 时它等于设计矩阵的数值秩；随着 \(\lambda\) 增大，各方向贡献逐渐减小。有效自由度描述模型的拟合灵活度，不是非零系数个数。

14.3.2 稳定实现与计算诊断

ridge_fit <- function(X, y, lambda, singular_tol = 1e-10) {
  if (length(lambda) != 1L || !is.finite(lambda) || lambda < 0) {
    stop("lambda 必须是非负有限数。")
  }

  prep <- standardize_xy(X, y)
  n <- nrow(prep$X)
  sv <- svd(
    prep$X,
    nu = min(dim(prep$X)),
    nv = min(dim(prep$X))
  )

  if (lambda == 0) {
    keep <- sv$d > singular_tol * max(sv$d)
    multiplier <- numeric(length(sv$d))
    multiplier[keep] <- 1 / sv$d[keep]
    effective_df <- sum(keep)
  } else {
    multiplier <- sv$d / (sv$d^2 + n * lambda)
    effective_df <- sum(sv$d^2 / (sv$d^2 + n * lambda))
  }

  beta_standardized <- drop(
    sv$v %*% (multiplier * crossprod(sv$u, prep$y))
  )
  restored <- restore_coefficients(beta_standardized, prep)
  residual_standardized <- prep$y - drop(prep$X %*% beta_standardized)

  list(
    intercept = restored$intercept,
    beta = restored$beta,
    beta_standardized = beta_standardized,
    lambda = lambda,
    effective_df = effective_df,
    objective = 0.5 * mean(residual_standardized^2) +
      0.5 * lambda * sum(beta_standardized^2),
    rank = sum(sv$d > singular_tol * max(sv$d)),
    prep = prep
  )
}

lambda = 0 时函数返回 SVD 广义逆对应的最小范数最小二乘解；lambda > 0 时使用 Ridge 收缩因子。函数同时返回有效自由度和目标函数值，便于核对不同惩罚强度下的计算。

ridge_demo_lambda <- c(0, 0.01, 0.1, 1, 10)
ridge_demo <- lapply(
  ridge_demo_lambda,
  function(lambda) ridge_fit(X_train, y_train, lambda)
)

ridge_demo_table <- data.frame(
  lambda = ridge_demo_lambda,
  effective_df = vapply(
    ridge_demo, function(fit) fit$effective_df, numeric(1)
  ),
  coefficient_norm = vapply(
    ridge_demo,
    function(fit) sqrt(sum(fit$beta_standardized^2)),
    numeric(1)
  ),
  training_RMSE = vapply(
    ridge_demo,
    function(fit) {
      sqrt(mean((y_train - predict_penalized(fit, X_train))^2))
    },
    numeric(1)
  )
)

knitr::kable(
  ridge_demo_table,
  row.names = FALSE,
  digits = 4,
  col.names = c(
    "$\\lambda$", "有效自由度", "标准化系数范数", "训练 RMSE"
  ),
  escape = FALSE
)
\(\lambda\) 有效自由度 标准化系数范数 训练 RMSE
0.00 60.000 2.7378 1.329
0.01 58.709 2.6620 1.330
0.10 50.209 2.1935 1.371
1.00 24.605 1.0439 1.792
10.00 4.947 0.2206 2.402

惩罚增强时，系数范数和有效自由度下降，训练误差则上升。训练误差上升本身不是缺点；是否值得接受这部分拟合损失，要由交叉验证和独立测试表现判断。

14.4 Lasso 与 Elastic Net：稀疏解的计算

14.4.1 从一范数惩罚到软阈值

令式 (14.1) 中 \(\alpha=1\)，Lasso 估计为 (Tibshirani 1996)

\[ \widehat{\boldsymbol\beta}_{\mathrm{lasso}}(\lambda) =\arg\min_{\boldsymbol\beta} \left\{ \frac{1}{2n}\lVert\mathbf y_{\mathrm{c}}-\mathbf X_{\mathrm{s}}\boldsymbol\beta\rVert_2^2 +\lambda\lVert\boldsymbol\beta\rVert_1 \right\}. \]

二范数平方在 0 点光滑，Ridge 系数通常只会连续接近 0。一范数在 0 点不可导，其最优性条件允许一整段梯度值对应 \(\beta_j=0\)，因而 Lasso 可以得到精确的零系数。

先考虑一维问题

\[ \min_b\left\{\frac{1}{2}(b-z)^2+\lambda|b|\right\}. \]

\(|b|\) 在 \(b=0\) 处的次梯度为区间 \([-1,1]\)。最优性条件

\[ 0\in b-z+\lambda\,\partial|b| \]

给出软阈值解

\[ S(z,\lambda) =\operatorname{sign}(z)(|z|-\lambda)_+, \qquad (a)_+=\max(a,0). \]

当 \(|z|\leq\lambda\) 时，结果为 0；超过阈值后，结果还会向 0 收缩 \(\lambda\)。硬阈值只删除小值而保留大值不变，软阈值则同时完成筛选和收缩。它也是一范数的近端算子，第 6.8 节已经介绍了这一观点，进一步讨论可参见 Parikh and Boyd (2014)。

soft_threshold <- function(z, threshold) {
  if (any(!is.finite(threshold)) || any(threshold < 0)) {
    stop("threshold 必须是非负有限数。")
  }
  sign(z) * pmax(abs(z) - threshold, 0)
}

soft_threshold(
  z = c(-2, -0.5, 0, 0.5, 2),
  threshold = 0.75
)
#> [1] -1.25  0.00  0.00  0.00  1.25

14.4.2 坐标下降与 KKT 诊断

坐标下降每次固定其他系数，只更新一个 \(\beta_j\)。记不含第 \(j\) 个变量贡献的部分残差为

\[ \mathbf r_j =\mathbf y_{\mathrm{c}}- \sum_{\substack{\ell=1\\ \ell\neq j}}^p \mathbf x_{\mathrm{s},\ell}\beta_\ell, \]

并定义

\[ z_j=\frac{1}{n}\mathbf x_{\mathrm{s},j}^\top\mathbf r_j, \qquad a_j=\frac{1}{n}\mathbf x_{\mathrm{s},j}^\top\mathbf x_{\mathrm{s},j}. \]

Lasso 的坐标更新为

\[ \beta_j\leftarrow\frac{S(z_j,\lambda)}{a_j}. \]

按本章标准化约定，\(a_j=1\)；代码仍显式计算它，以便暴露预处理错误。高效实现不必在每次更新时重新计算完整矩阵乘法。若当前残差为\(\mathbf r=\mathbf y_{\mathrm{c}}-\mathbf X_{\mathrm{s}}\boldsymbol\beta\)，则先加回旧的第 \(j\) 项得到\(\mathbf r_j=\mathbf r+\mathbf x_{\mathrm{s},j}\beta_j\)，更新系数后再令\(\mathbf r=\mathbf r_j-\mathbf x_{\mathrm{s},j}\beta_j^{\mathrm{new}}\)。一次更新只访问一列。

坐标更新给出了迭代方法，但还没有说明何时可以停止。“相邻两次系数变化很小”并不足以证明已经接近最优解。令光滑损失的梯度为

\[ g_j(\boldsymbol\beta) =\frac{1}{n}\mathbf x_{\mathrm{s},j}^\top (\mathbf X_{\mathrm{s}}\boldsymbol\beta-\mathbf y_{\mathrm{c}}). \]

Lasso 的 Karush–Kuhn–Tucker（KKT）条件为

\[ \begin{cases} g_j(\widehat{\boldsymbol\beta})+\lambda\operatorname{sign}(\hat\beta_j)=0, &\hat\beta_j\neq0,\\ |g_j(\widehat{\boldsymbol\beta})|\leq\lambda, &\hat\beta_j=0. \end{cases} \]

因此可以把所有坐标违反上述条件的最大程度作为计算诊断。下面的函数维护残差，并同时检查相对系数变化和 KKT 违反量。

elastic_net_cd <- function(X, y, lambda, alpha = 1,
                           beta_init = NULL, max_iter = 5000L,
                           coefficient_tol = 1e-8,
                           kkt_tol = 1e-7) {
  if (length(lambda) != 1L || !is.finite(lambda) || lambda < 0) {
    stop("lambda 必须是非负有限数。")
  }
  if (length(alpha) != 1L || !is.finite(alpha) ||
      alpha < 0 || alpha > 1) {
    stop("alpha 必须位于 0 与 1 之间。")
  }

  prep <- standardize_xy(X, y)
  Xs <- prep$X
  yc <- prep$y
  n <- nrow(Xs)
  p <- ncol(Xs)
  x2 <- colMeans(Xs^2)

  if (is.null(beta_init)) {
    beta <- numeric(p)
  } else {
    beta <- as.numeric(beta_init)
    if (length(beta) != p || any(!is.finite(beta))) {
      stop("beta_init 必须是长度等于 X 列数的有限向量。")
    }
  }

  residual <- yc - drop(Xs %*% beta)
  objective_trace <- numeric(max_iter)
  kkt_trace <- numeric(max_iter)
  converged <- FALSE

  for (iter in seq_len(max_iter)) {
    beta_old <- beta

    for (j in seq_len(p)) {
      partial_residual <- residual + Xs[, j] * beta[j]
      z_j <- mean(Xs[, j] * partial_residual)
      beta_new <- soft_threshold(z_j, lambda * alpha) /
        (x2[j] + lambda * (1 - alpha))
      residual <- partial_residual - Xs[, j] * beta_new
      beta[j] <- beta_new
    }

    gradient <- -drop(crossprod(Xs, residual)) / n
    active <- abs(beta) > 1e-12
    violation <- numeric(p)
    violation[active] <- abs(
      gradient[active] + lambda * (1 - alpha) * beta[active] +
        lambda * alpha * sign(beta[active])
    )
    violation[!active] <- pmax(
      abs(gradient[!active]) - lambda * alpha,
      0
    )

    kkt_violation <- max(violation)
    relative_change <- max(abs(beta - beta_old)) /
      (1 + max(abs(beta_old)))
    objective_trace[iter] <- 0.5 * mean(residual^2) +
      lambda * (
        alpha * sum(abs(beta)) +
          0.5 * (1 - alpha) * sum(beta^2)
      )
    kkt_trace[iter] <- kkt_violation

    if (relative_change <= coefficient_tol &&
        kkt_violation <= kkt_tol) {
      converged <- TRUE
      break
    }
  }

  restored <- restore_coefficients(beta, prep)
  list(
    intercept = restored$intercept,
    beta = restored$beta,
    beta_standardized = beta,
    lambda = lambda,
    alpha = alpha,
    converged = converged,
    iterations = iter,
    kkt_violation = kkt_violation,
    objective = objective_trace[iter],
    objective_trace = objective_trace[seq_len(iter)],
    kkt_trace = kkt_trace[seq_len(iter)],
    prep = prep
  )
}

lasso_fit <- function(X, y, lambda, ...) {
  elastic_net_cd(X, y, lambda, alpha = 1, ...)
}

如果函数达到 max_iter 却未满足条件，converged 会返回 FALSE。此时不应只增加输出位数或忽略警告，而应检查迭代上限、容差、变量尺度以及问题是否存在高度相关或非唯一解。

14.4.3 正则化路径与暖启动

Lasso 的零向量满足 KKT 条件的最小惩罚强度为

\[ \lambda_{\max} =\max_{1\leq j\leq p}\left| \frac{1}{n}\mathbf x_{\mathrm{s},j}^\top\mathbf y_{\mathrm{c}} \right|. \]

从 \(\lambda_{\max}\) 开始逐渐减小 \(\lambda\)，可以得到从全零系数到较复杂模型的正则化路径。相邻\(\lambda\) 的解通常接近，因此把前一个解作为下一个问题的初值，即暖启动（warm start），可以减少迭代次数。

lasso_lambda_max <- function(X, y) {
  prep <- standardize_xy(X, y)
  max(abs(drop(crossprod(prep$X, prep$y))) / nrow(prep$X))
}

lambda_max_full <- lasso_lambda_max(X_train, y_train)
lambda_fraction <- exp(seq(log(1), log(0.01), length.out = 40))
lambda_path <- lambda_max_full * lambda_fraction

lasso_path_beta <- matrix(
  NA_real_,
  nrow = length(lambda_path),
  ncol = ncol(X_train)
)
lasso_path_iterations <- integer(length(lambda_path))
beta_start <- NULL

for (i in seq_along(lambda_path)) {
  fit_i <- lasso_fit(
    X_train,
    y_train,
    lambda = lambda_path[i],
    beta_init = beta_start
  )
  stopifnot(fit_i$converged)
  lasso_path_beta[i, ] <- fit_i$beta_standardized
  lasso_path_iterations[i] <- fit_i$iterations
  beta_start <- fit_i$beta_standardized
}

matplot(
  log10(lambda_path),
  lasso_path_beta,
  type = "l",
  lty = 1,
  col = c(rep("firebrick", 5), rep("grey75", ncol(X_train) - 5)),
  xlab = expression(log[10](lambda)),
  ylab = "标准化系数"
)
abline(h = 0, col = "grey40", lty = 3)
Lasso标准化系数随惩罚强度变化的路径

图14.1: Lasso标准化系数随惩罚强度变化的路径

图中红线对应具有真实信号的前 5 个变量，灰线对应无关变量。路径并不保证先进入的变量就是真实变量，尤其在强相关变量之间，样本扰动可能改变进入顺序。正则化路径展示的是给定数据上的优化结果，变量选择的稳定性仍需用重复抽样或外部数据检查。

14.4.4 Elastic Net：在稀疏与稳定之间

当解释变量高度相关时，Lasso 可能从一组相似变量中任意保留少数变量，选择结果随样本变化。ElasticNet 在一范数惩罚之外加入二次惩罚 (Zou and Hastie 2005)。其目标正是式(14.1)；坐标更新为

\[ \beta_j\leftarrow \frac{S(z_j,\lambda\alpha)} {a_j+\lambda(1-\alpha)}. \]

分子中的软阈值产生稀疏性，分母中的二次项稳定相关方向。elastic_net_cd() 已经按照这一公式实现：alpha = 1 为 Lasso，alpha = 0 为坐标下降版本的 Ridge，中间值为 Elastic Net。实践中 \(\alpha\)也是调节参数；若同时搜索 \(\alpha\) 和 \(\lambda\)，两层选择都必须包含在交叉验证过程内。

Ridge、Lasso 与 Elastic Net 的主要差异可概括如下。

性质 Ridge Lasso Elastic Net
系数是否可为精确 0 通常不会 可以 可以
相关变量的处理 倾向共同收缩 可能只保留其中少数 较容易成组保留并收缩
解的唯一性 \(\lambda>0\) 时唯一 可能因设计矩阵而不唯一 含正二次惩罚时通常更稳定
主要计算 SVD、线性方程组 软阈值、坐标下降 带二次项的坐标下降

稀疏系数便于存储和预测，但“系数等于 0”是惩罚与样本共同产生的结果，不等同于变量在总体中完全无效，更不能自动获得因果解释。

14.5 ADMM：另一种 Lasso 分解方法（选读）

坐标下降利用目标函数对单个坐标的结构。交替方向乘子法（alternating direction method ofmultipliers，ADMM）则通过变量分裂，把光滑拟合项与不可导惩罚项交给不同变量处理。令\(\boldsymbol\beta=\mathbf z\)，Lasso 可写为

\[ \begin{aligned} \min_{\boldsymbol\beta,\mathbf z}\quad &\frac{1}{2n} \lVert\mathbf y_{\mathrm{c}}-\mathbf X_{\mathrm{s}}\boldsymbol\beta\rVert_2^2 +\lambda\lVert\mathbf z\rVert_1,\\ \text{subject to}\quad&\boldsymbol\beta=\mathbf z. \end{aligned} \]

使用缩放对偶变量 \(\mathbf u\) 和算法参数 \(\rho>0\)，三步更新为

\[ \begin{aligned} \boldsymbol\beta^{(t+1)} &=\left(\frac{\mathbf X_{\mathrm{s}}^\top\mathbf X_{\mathrm{s}}}{n}+\rho\mathbf I_p\right)^{-1} \left\{\frac{\mathbf X_{\mathrm{s}}^\top\mathbf y_{\mathrm{c}}}{n} +\rho(\mathbf z^{(t)}-\mathbf u^{(t)})\right\},\\ \mathbf z^{(t+1)} &=S\left(\boldsymbol\beta^{(t+1)}+\mathbf u^{(t)}, \frac{\lambda}{\rho}\right),\\ \mathbf u^{(t+1)} &=\mathbf u^{(t)}+\boldsymbol\beta^{(t+1)}-\mathbf z^{(t+1)}. \end{aligned} \]

第一步是带二次稳定项的线性方程组，第二步把软阈值逐分量作用于向量，第三步更新约束的累计偏差。\(\rho\) 影响迭代路径和收敛速度，却不是模型的正则化强度。关于 ADMM 的系统推导可参见 Boyd et al. (2011)。

停止条件同时检查原始残差和对偶残差：

\[ \mathbf r^{(t)}=\boldsymbol\beta^{(t)}-\mathbf z^{(t)}, \qquad \mathbf s^{(t)}=\rho(\mathbf z^{(t)}-\mathbf z^{(t-1)}). \]

固定的绝对阈值会随参数维数和系数尺度改变含义。下面按照绝对容差与相对容差构造维数相关的停止阈值，并缓存第一步所需矩阵的 Cholesky 分解。

lasso_admm <- function(X, y, lambda, rho = 1,
                       max_iter = 5000L,
                       abs_tol = 1e-6,
                       rel_tol = 1e-5) {
  if (length(lambda) != 1L || !is.finite(lambda) || lambda < 0) {
    stop("lambda 必须是非负有限数。")
  }
  if (length(rho) != 1L || !is.finite(rho) || rho <= 0) {
    stop("rho 必须是正的有限数。")
  }

  prep <- standardize_xy(X, y)
  Xs <- prep$X
  yc <- prep$y
  n <- nrow(Xs)
  p <- ncol(Xs)

  beta <- numeric(p)
  z <- numeric(p)
  u <- numeric(p)
  system_matrix <- crossprod(Xs) / n + rho * diag(p)
  chol_factor <- chol(system_matrix)
  xy <- drop(crossprod(Xs, yc)) / n

  solve_system <- function(rhs) {
    backsolve(
      chol_factor,
      forwardsolve(t(chol_factor), rhs)
    )
  }

  converged <- FALSE
  for (iter in seq_len(max_iter)) {
    z_old <- z
    beta <- drop(solve_system(xy + rho * (z - u)))
    z <- soft_threshold(beta + u, lambda / rho)
    u <- u + beta - z

    primal_residual <- sqrt(sum((beta - z)^2))
    dual_residual <- rho * sqrt(sum((z - z_old)^2))
    primal_tolerance <- sqrt(p) * abs_tol +
      rel_tol * max(sqrt(sum(beta^2)), sqrt(sum(z^2)))
    dual_tolerance <- sqrt(p) * abs_tol +
      rel_tol * rho * sqrt(sum(u^2))

    if (primal_residual <= primal_tolerance &&
        dual_residual <= dual_tolerance) {
      converged <- TRUE
      break
    }
  }

  restored <- restore_coefficients(z, prep)
  fitted_residual <- yc - drop(Xs %*% z)
  list(
    intercept = restored$intercept,
    beta = restored$beta,
    beta_standardized = z,
    lambda = lambda,
    rho = rho,
    converged = converged,
    iterations = iter,
    primal_residual = primal_residual,
    dual_residual = dual_residual,
    primal_tolerance = primal_tolerance,
    dual_tolerance = dual_tolerance,
    objective = 0.5 * mean(fitted_residual^2) +
      lambda * sum(abs(z)),
    prep = prep
  )
}

下面在同一个 \(\lambda\) 上比较坐标下降和 ADMM。凸目标的最优解应当一致；“两段代码都运行结束”不是核对，系数差、目标函数差和各自的停止诊断才是证据。

lambda_check <- 0.15 * lambda_max_full
fit_cd_check <- lasso_fit(
  X_train, y_train, lambda = lambda_check,
  kkt_tol = 1e-8
)
fit_admm_check <- lasso_admm(
  X_train, y_train, lambda = lambda_check,
  rho = 1,
  abs_tol = 1e-7,
  rel_tol = 1e-6
)

stopifnot(
  fit_cd_check$converged,
  fit_admm_check$converged,
  max(abs(
    fit_cd_check$beta_standardized -
      fit_admm_check$beta_standardized
  )) < 1e-4
)

admm_check <- data.frame(
  quantity = c(
    "坐标下降是否收敛",
    "ADMM是否收敛",
    "标准化系数最大绝对差",
    "目标函数绝对差",
    "坐标下降KKT违反量",
    "ADMM原始残差",
    "ADMM对偶残差"
  ),
  result = c(
    as.character(fit_cd_check$converged),
    as.character(fit_admm_check$converged),
    formatC(max(abs(
      fit_cd_check$beta_standardized -
        fit_admm_check$beta_standardized
    )), format = "e", digits = 3),
    formatC(
      abs(fit_cd_check$objective - fit_admm_check$objective),
      format = "e", digits = 3
    ),
    formatC(fit_cd_check$kkt_violation, format = "e", digits = 3),
    formatC(fit_admm_check$primal_residual, format = "e", digits = 3),
    formatC(fit_admm_check$dual_residual, format = "e", digits = 3)
  )
)

knitr::kable(
  admm_check,
  row.names = FALSE,
  col.names = c("核对量", "结果")
)
核对量 结果
坐标下降是否收敛 TRUE
ADMM是否收敛 TRUE
标准化系数最大绝对差 3.764e-06
目标函数绝对差 7.632e-12
坐标下降KKT违反量 3.456e-09
ADMM原始残差 1.511e-06
ADMM对偶残差 1.755e-06

对于普通 Lasso，成熟坐标下降实现通常很高效。ADMM 的价值更多体现在带额外约束、分组结构或可分解目标的问题。当前实现显式形成 \(\mathbf X_{\mathrm{s}}^\top\mathbf X_{\mathrm{s}}\)，若 \(p\) 很大或设计矩阵稀疏，这一步可能成为内存瓶颈；变量分裂并不自动保证每个线性代数步骤都适合大规模计算。

本章函数用于展示目标函数、更新式和停止条件之间的联系，不应被理解为生产级求解器。实际分析宜使用经过系统测试、能够处理稀疏矩阵并实现筛选规则的成熟软件，例如 glmnet；但即使调用成熟软件，仍要核对目标函数缩放、标准化约定和收敛信息，不能只比较输出中的 \(\lambda\) 数值。

14.6 调节参数选择与最终评价

14.6.1 折内预处理与惩罚网格

\(\lambda\) 控制模型从弱惩罚到强惩罚的路径，\(\alpha\) 决定一范数与二次惩罚的比例。两者都属于模型选择的一部分。按照第 13 章的原则，每个验证折都要执行一条完整管道：

  1. 只用当前训练折估计中心、尺度和惩罚路径；
  2. 在训练折上拟合各候选参数；
  3. 使用训练折的中心和尺度对验证折产生预测；
  4. 保存每个观测的折外损失；
  5. 所有模型共用同一组分折，再比较汇总误差。

Lasso 在给定训练数据上的 \(\lambda_{\max}\) 依赖\(\mathbf X_{\mathrm{s}}^\top\mathbf y_{\mathrm{c}}\)。如果先用全样本响应计算\(\lambda_{\max}\)，再进行交叉验证，验证折已经影响了候选路径。下面对 Lasso 和 Elastic Net 使用\(\lambda/\lambda_{\max}\) 的固定比例网格，并在每个训练折内部重新计算 \(\lambda_{\max}\)。Ridge 没有使所有系数恰好为 0 的有限 \(\lambda_{\max}\)，因此预先给出一个从强惩罚到弱惩罚的固定网格。按本章的目标函数和自变量标准化约定，Ridge 的 \(\lambda\) 不随响应变量单位的线性变换而改变；Lasso 的\(\lambda_{\max}\) 则会随响应尺度同比例变化。Elastic Net 同时混合一次和二次齐次的惩罚，改变响应尺度还可能改变两部分的相对作用，因此比较结果时必须固定预处理与软件约定。这也是必须说明目标函数尺度的一个具体原因。

cv_penalized <- function(X, y, fold_id, alpha, penalty_grid) {
  X <- as.matrix(X)
  y <- as.numeric(y)
  n <- nrow(X)
  fold_levels <- sort(unique(fold_id))

  if (length(y) != n || length(fold_id) != n || anyNA(fold_id)) {
    stop("X、y 与 fold_id 的长度不一致。")
  }
  if (length(fold_levels) < 2L) {
    stop("交叉验证至少需要两个分折。")
  }
  if (length(alpha) != 1L || !is.finite(alpha) ||
      alpha < 0 || alpha > 1) {
    stop("alpha 必须位于 0 与 1 之间。")
  }
  if (any(!is.finite(penalty_grid)) || any(penalty_grid <= 0) ||
      is.unsorted(-penalty_grid)) {
    stop("penalty_grid 必须按从强到弱的顺序递减排列。")
  }

  n_grid <- length(penalty_grid)
  squared_error <- matrix(NA_real_, nrow = n, ncol = n_grid)
  convergence <- matrix(
    TRUE,
    nrow = length(fold_levels),
    ncol = n_grid
  )

  for (fold_position in seq_along(fold_levels)) {
    in_validation <- fold_id == fold_levels[fold_position]
    X_fit <- X[!in_validation, , drop = FALSE]
    y_fit <- y[!in_validation]
    X_validation <- X[in_validation, , drop = FALSE]

    if (alpha == 0) {
      lambda_values <- penalty_grid
    } else {
      lambda_values <- penalty_grid *
        lasso_lambda_max(X_fit, y_fit) / alpha
    }

    beta_start <- NULL
    for (j in seq_len(n_grid)) {
      if (alpha == 0) {
        fit <- ridge_fit(X_fit, y_fit, lambda_values[j])
      } else {
        fit <- elastic_net_cd(
          X_fit,
          y_fit,
          lambda = lambda_values[j],
          alpha = alpha,
          beta_init = beta_start
        )
        beta_start <- fit$beta_standardized
        convergence[fold_position, j] <- fit$converged
      }

      prediction <- predict_penalized(fit, X_validation)
      squared_error[in_validation, j] <-
        (y[in_validation] - prediction)^2
    }
  }

  if (any(!is.finite(squared_error))) {
    stop("交叉验证损失包含非有限值。")
  }
  if (!all(convergence)) {
    stop("至少一个训练折上的坐标下降未收敛。")
  }

  fold_mse <- vapply(
    seq_len(n_grid),
    function(j) tapply(
      squared_error[, j], fold_id, mean
    ),
    numeric(length(fold_levels))
  )
  mean_mse <- colMeans(squared_error)
  fold_se <- apply(fold_mse, 2, sd) / sqrt(length(fold_levels))
  index_min <- which.min(mean_mse)
  one_se_limit <- mean_mse[index_min] + fold_se[index_min]
  index_1se <- which(mean_mse <= one_se_limit)[1]

  list(
    alpha = alpha,
    penalty_grid = penalty_grid,
    mean_mse = mean_mse,
    fold_se = fold_se,
    index_min = index_min,
    index_1se = index_1se,
    one_se_limit = one_se_limit,
    squared_error = squared_error,
    fold_mse = fold_mse
  )
}

penalty_grid 按惩罚从强到弱排列，所以满足一标准误阈值的第一个候选值就是较强惩罚。这里的折间标准误是常见的一标准误规则所用近似；各训练集高度重叠，不能把它当作严格的独立重复抽样标准误。一标准误规则的作用是：在验证误差与最小值难以区分时，倾向选择更简单、惩罚更强的模型。

14.6.2 用共同分折比较三类惩罚

下面使用固定的五折划分。Ridge 的网格是绝对 \(\lambda\)，Lasso 与 Elastic Net 的网格是\(\lambda/\lambda_{\max}\)；这些数值不能直接横向比较，能够比较的是使用同一折外观测计算的预测误差。

set.seed(1401)
K <- 5L
fold_id <- sample(rep(seq_len(K), length.out = nrow(X_train)))

ridge_grid <- 10^seq(2, -3, length.out = 36)
relative_grid <- exp(seq(log(1), log(0.01), length.out = 36))

cv_ridge <- cv_penalized(
  X_train, y_train, fold_id,
  alpha = 0,
  penalty_grid = ridge_grid
)
cv_elastic <- cv_penalized(
  X_train, y_train, fold_id,
  alpha = 0.5,
  penalty_grid = relative_grid
)
cv_lasso <- cv_penalized(
  X_train, y_train, fold_id,
  alpha = 1,
  penalty_grid = relative_grid
)

cv_objects <- list(
  Ridge = cv_ridge,
  `Elastic Net` = cv_elastic,
  Lasso = cv_lasso
)

old_par <- par(mfrow = c(1, 3), mar = c(4.2, 4, 2, 0.8))
for (method in names(cv_objects)) {
  result <- cv_objects[[method]]
  plot(
    log10(result$penalty_grid),
    sqrt(result$mean_mse),
    type = "l",
    lwd = 2,
    xlab = if (result$alpha == 0) {
      expression(log[10](lambda))
    } else {
      expression(log[10](lambda / lambda[max]))
    },
    ylab = "CV RMSE",
    main = method
  )
  points(
    log10(result$penalty_grid[result$index_min]),
    sqrt(result$mean_mse[result$index_min]),
    pch = 19,
    col = "steelblue"
  )
  points(
    log10(result$penalty_grid[result$index_1se]),
    sqrt(result$mean_mse[result$index_1se]),
    pch = 17,
    col = "firebrick"
  )
}
Ridge、Elastic Net与Lasso的五折交叉验证曲线

图14.2: Ridge、Elastic Net与Lasso的五折交叉验证曲线 

par(old_par)

蓝色圆点表示最小交叉验证误差，红色三角表示一标准误规则选择的较强惩罚。候选网格必须覆盖曲线最低点附近；如果最优值落在网格边界，应扩展网格后重新计算，而不是直接接受边界结果。

14.6.3 全训练数据重拟合与测试评价

交叉验证选定参数后，要在全部训练数据上重新拟合。测试集此前没有参与标准化、路径构造、参数选择或一标准误规则。下面采用一标准误选择，并在测试集上评价一次。

fit_selected_penalty <- function(X, y, cv_result) {
  selected_grid_value <-
    cv_result$penalty_grid[cv_result$index_1se]

  if (cv_result$alpha == 0) {
    selected_lambda <- selected_grid_value
    fit <- ridge_fit(X, y, selected_lambda)
  } else {
    selected_lambda <- selected_grid_value *
      lasso_lambda_max(X, y) / cv_result$alpha
    fit <- elastic_net_cd(
      X,
      y,
      lambda = selected_lambda,
      alpha = cv_result$alpha
    )
  }

  if (!is.null(fit$converged) && !fit$converged) {
    stop("全训练数据上的最终拟合未收敛。")
  }
  fit$selected_grid_value <- selected_grid_value
  fit
}

selected_fits <- Map(
  function(result) fit_selected_penalty(X_train, y_train, result),
  cv_objects
)

ols_fit <- ridge_fit(X_train, y_train, lambda = 0)
all_final_fits <- c(list(OLS = ols_fit), selected_fits)

test_comparison <- do.call(
  rbind,
  lapply(names(all_final_fits), function(method) {
    fit <- all_final_fits[[method]]
    train_prediction <- predict_penalized(fit, X_train)
    test_prediction <- predict_penalized(fit, X_test)
    data.frame(
      method = method,
      lambda = if (method == "OLS") 0 else fit$lambda,
      nonzero = sum(abs(fit$beta) > 1e-8),
      training_RMSE = sqrt(mean((y_train - train_prediction)^2)),
      test_RMSE = sqrt(mean((y_test - test_prediction)^2)),
      test_MAE = mean(abs(y_test - test_prediction))
    )
  })
)

knitr::kable(
  test_comparison,
  row.names = FALSE,
  digits = 4,
  col.names = c(
    "方法", "$\\lambda$", "非零系数", "训练 RMSE",
    "测试 RMSE", "测试 MAE"
  ),
  escape = FALSE
)
方法 \(\lambda\) 非零系数 训练 RMSE 测试 RMSE 测试 MAE
OLS 0.0000 60 1.329 2.030 1.618
Ridge 1.3895 60 1.894 2.172 1.748
Elastic Net 0.3448 19 1.727 1.786 1.435
Lasso 0.2558 9 1.776 1.714 1.375

OLS 的训练误差通常最小，却需要估计全部系数；正则化模型接受一定的训练拟合损失，以换取较稳定的样本外预测。Lasso 和 Elastic Net 的非零系数数目反映稀疏程度，Ridge 的系数通常全部非零，复杂度应结合有效自由度理解。一次模拟结果不能证明某种方法总是更优；变量相关结构、信号稀疏程度、噪声大小和样本量都会改变比较结果。

如果还要从多个 \(\alpha\) 中选择 Elastic Net，应把所有候选 \(\alpha\) 放进同一交叉验证。若没有独立测试集而又要报告整个调参流程的性能，应使用第 13 章介绍的嵌套交叉验证。

14.7 从正则化解到统计结论

14.7.1 系数解释与选择稳定性

正则化模型的系数表需要结合惩罚与预处理解释。

第一，原始尺度系数和标准化尺度系数回答不同问题。原始尺度系数适合产生预测和说明变量单位变化的影响；标准化系数便于比较惩罚路径，但其大小仍受变量分布和相关结构影响，不能简单排序为“重要性”。

第二，Lasso 选择的是一组有利于当前预测目标的变量。在相关变量中，多个不同变量集合可能产生近似相同的预测；轻微改变样本或 \(\lambda\) 就可能改变非零集合。Elastic Net 往往比纯 Lasso 更容易共同保留相关变量，但也不能保证选择结果稳定。

第三，选择后的常规推断需要特别处理。先用同一数据选择变量，再把最终非零变量当作事先指定，直接套用普通最小二乘标准误和 \(p\) 值，会忽略选择步骤的不确定性。若目标是确认变量效应，应依靠研究设计、独立样本、样本拆分或专门的选择后推断方法，而不能把预测型正则化自动解释为因果发现。

第四，数值上的“非零”需要容差。迭代算法可能返回 \(10^{-12}\) 量级的系数；报告非零个数时应说明阈值，并同时检查 KKT 违反量。只把很小的数打印成 0，不代表最优性条件已经满足。

下面固定原数据选定的 Lasso 路径比例，在 Bootstrap 样本内重新计算 \(\lambda_{\max}\) 并拟合，用变量被选中的频率描述样本扰动下的选择稳定性。它是诊断量，不是变量具有真实效应的后验概率或显著性概率。

set.seed(1402)
B <- 50L
lasso_fraction_1se <-
  cv_lasso$penalty_grid[cv_lasso$index_1se]
selection_record <- matrix(
  FALSE,
  nrow = B,
  ncol = ncol(X_train),
  dimnames = list(NULL, colnames(X_train))
)

for (b in seq_len(B)) {
  bootstrap_id <- sample(
    seq_len(nrow(X_train)),
    replace = TRUE
  )
  X_boot <- X_train[bootstrap_id, , drop = FALSE]
  y_boot <- y_train[bootstrap_id]
  lambda_boot <- lasso_fraction_1se *
    lasso_lambda_max(X_boot, y_boot)
  fit_boot <- lasso_fit(X_boot, y_boot, lambda_boot)
  stopifnot(fit_boot$converged)
  selection_record[b, ] <- abs(fit_boot$beta) > 1e-8
}

selection_frequency <- colMeans(selection_record)
top_selection <- order(selection_frequency, decreasing = TRUE)[1:12]
stability_table <- data.frame(
  variable = names(selection_frequency)[top_selection],
  true_coefficient = beta_true[top_selection],
  selection_frequency = selection_frequency[top_selection]
)

knitr::kable(
  stability_table,
  row.names = FALSE,
  digits = 3,
  col.names = c("变量", "真实系数", "Bootstrap选择频率")
)
变量 真实系数 Bootstrap选择频率
x1 1.2 1.00
x2 -1.0 1.00
x4 0.6 1.00
x3 0.8 0.96
x5 -0.5 0.96
x41 0.0 0.72
x57 0.0 0.72
x51 0.0 0.62
x14 0.0 0.50
x47 0.0 0.44
x34 0.0 0.42
x38 0.0 0.42

真实信号较强的变量通常有较高选择频率，但相关的无关变量也可能被选中，较弱真实变量也可能在部分样本中消失。这个现象说明“预测表现稳定”和“变量集合稳定”是两个不同要求。

14.7.2 可核验的计算证据

正则化分析的统计结论依赖数据划分、预处理、参数选择和数值优化。为了使结果可以核验，报告应覆盖整个计算链条。

计算环节 应保留的证据
数据输入 训练、验证和测试数据的定义；时间、分组或空间边界；模型矩阵与缺失值处理
预处理与目标 中心和尺度；目标函数的缩放；截距是否惩罚；\(\lambda\) 与 \(\alpha\) 的定义
参数选择 候选网格、分折编号、选择规则和随机种子
数值求解 收敛状态、迭代次数、KKT 违反量或原始—对偶残差；软件版本与矩阵表示
统计结果 非零判定阈值、相关变量结构与选择稳定性；折外或测试指标及简单基准模型

只报告最终系数和一个 \(\lambda\)，无法判断预处理是否泄漏、调参是否公平，也无法核对优化是否收敛。

14.8 稀疏计算

14.8.1 统计稀疏性与存储稀疏性

高维正则化中常见两种稀疏性，二者不能混为一谈。

类型 含义 主要计算收益
系数稀疏 \(\widehat{\boldsymbol\beta}\) 中许多分量为 0 减少预测所需变量，便于存储模型
设计矩阵稀疏 \(\mathbf X\) 中大多数元素为 0 减少数据存储和矩阵—向量运算成本

Lasso 可能产生稀疏系数，但输入的 \(\mathbf X\) 不一定稀疏；文本词频、用户—商品记录和大量分类变量生成的哑变量设计矩阵则常常稀疏，即使所用模型不产生零系数。算法只有在数据结构和运算过程中保持这种结构，才能得到实际的内存和速度收益。

若 \(n\times p\) 矩阵只有 \(m\) 个非零元素，其密度为

\[ \operatorname{density}(\mathbf X)=\frac{m}{np}. \]

稠密表示保存全部 \(np\) 个位置；压缩稀疏列格式主要保存非零值、行索引和各列起点。当 \(m\ll np\) 时，两者的存储差别很大。下面构造一个稀疏购买记录矩阵。

set.seed(1403)
n_user <- 1000L
n_item <- 2000L
n_record <- 10000L

purchase_pair <- unique(data.frame(
  user = sample(n_user, n_record, replace = TRUE),
  item = sample(n_item, n_record, replace = TRUE)
))

purchase_sparse <- Matrix::sparseMatrix(
  i = purchase_pair$user,
  j = purchase_pair$item,
  x = 1,
  dims = c(n_user, n_item)
)
purchase_dense <- as.matrix(purchase_sparse)

sparse_density <- Matrix::nnzero(purchase_sparse) /
  prod(dim(purchase_sparse))
storage_table <- data.frame(
  representation = c("稠密矩阵", "稀疏矩阵"),
  size_MB = as.numeric(c(
    object.size(purchase_dense),
    object.size(purchase_sparse)
  )) / 1024^2,
  density = sparse_density
)

knitr::kable(
  storage_table,
  row.names = FALSE,
  digits = 4,
  col.names = c("表示方式", "对象大小（MB）", "非零元素比例")
)
表示方式 对象大小（MB） 非零元素比例
稠密矩阵 15.2590 0.005
稀疏矩阵 0.1232 0.005

对象大小还包含索引和类信息，因此不会严格等于简单的元素个数公式。关键结论是：低密度时，稀疏表示避免为大量 0 分配双精度存储。

14.8.2 稀疏表示如何改变计算路径

set.seed(1404)
item_score <- rnorm(n_item)
dense_result <- drop(purchase_dense %*% item_score)
sparse_result <- drop(purchase_sparse %*% item_score)

sparse_multiplication_check <- max(
  abs(dense_result - sparse_result)
)
stopifnot(sparse_multiplication_check < 1e-12)
sparse_multiplication_check
#> [1] 0

两种表示对应同一个线性变换，但稀疏乘法只需处理非零元素。坐标下降按列访问设计矩阵，与压缩列格式自然配合。相反，显式构造 \(\mathbf X^\top\mathbf X\) 可能产生填充：即使 \(\mathbf X\) 很稀疏，其交叉乘积也可能稠密。此时直接使用矩阵—向量乘法或迭代线性代数方法更合适。

中心化也可能破坏稀疏性。原来为 0 的元素减去非零列均值后都会变为非零。成熟软件通常单独记录中心和尺度，在运算中隐式完成校正，而不真正建立中心化后的稠密矩阵。本章的教学函数使用普通矩阵，便于展示算法；面对大规模稀疏输入时，应选用能够原生保留稀疏结构的软件实现。

分类变量经独热编码后也常产生稀疏设计矩阵。Matrix::sparse.model.matrix() 可以直接构造稀疏形式，避免先生成密集矩阵再转换。

set.seed(1405)
n_observation <- 2500L
industry <- factor(
  sample(paste0("industry_", 1:80), n_observation, replace = TRUE)
)
region <- factor(
  sample(paste0("region_", 1:50), n_observation, replace = TRUE)
)

X_dense_dummy <- model.matrix(~ industry + region)
X_sparse_dummy <- Matrix::sparse.model.matrix(~ industry + region)

dummy_table <- data.frame(
  representation = c("稠密设计矩阵", "稀疏设计矩阵"),
  rows = nrow(X_dense_dummy),
  columns = ncol(X_dense_dummy),
  size_MB = as.numeric(c(
    object.size(X_dense_dummy),
    object.size(X_sparse_dummy)
  )) / 1024^2,
  density = c(
    mean(X_dense_dummy != 0),
    Matrix::nnzero(X_sparse_dummy) / prod(dim(X_sparse_dummy))
  )
)

knitr::kable(
  dummy_table,
  row.names = FALSE,
  digits = 4,
  col.names = c(
    "表示方式", "行数", "列数", "对象大小（MB）", "非零元素比例"
  )
)
表示方式 行数 列数 对象大小（MB） 非零元素比例
稠密设计矩阵 2500 129 2.625 0.023
稀疏设计矩阵 2500 129 0.251 0.023

正则化还可以从概率模型的角度理解。后续章节将说明，高斯型先验与 Ridge 型惩罚、双指数型先验与Lasso 型惩罚之间存在对应关系。这里不提前使用尚未建立的后验公式，但需要注意：这种对应关系取决于似然和惩罚项的缩放；惩罚优化给出的是点估计，完整概率分析还要描述参数不确定性并把它传播到预测。第 15 章将从先验、似然和后验的基本概念开始讨论。

14.9 本章小结

正则化主动引入收缩偏差，以降低不稳定估计带来的方差；它是否改善预测取决于两者的合计变化。Ridge、Lasso 和 Elastic Net 可以写在统一的惩罚最小二乘框架下。截距通常不惩罚，自变量需要按照明确约定标准化；在交叉验证中，这些预处理量只能从训练折估计。目标函数的缩放决定 \(\lambda\) 的数值含义，跨软件比较时必须先核对尺度。

Ridge 通过缩小小奇异值方向的影响改善稳定性，SVD 揭示了收缩因子和有效自由度。Lasso 的零系数来自一范数的次梯度结构，软阈值是坐标下降的基本计算单元。维护残差可以减少重复运算，KKT 条件则提供比“系数几乎不动”更直接的最优性检查。Elastic Net 通过二次惩罚稳定相关变量；ADMM 从变量分裂的角度给出另一种计算路径。

调节参数属于模型选择，必须用折外预测评价。各方法应共享分折，Lasso 路径和标准化要在折内重算，独立测试集只在最终阶段使用。系数稀疏与设计矩阵稀疏是不同概念；后者只有在存储和运算都保持稀疏结构时才会转化为计算收益。完整报告还应保留参数网格、数据分折、收敛诊断和选择稳定性等计算证据。

14.10 思考题

  1. 式 (14.1) 若去掉 \(1/n\)，怎样换算 \(\lambda\) 才能保持同一个最优解？为什么不同软件输出的\(\lambda\) 不能只按数值直接比较？
  2. 推导 Ridge 的 SVD 表达式，并说明 \(d_j/(d_j^2+n\lambda)\) 如何抑制小奇异值方向。
  3. Ridge 的有效自由度为什么通常不是整数？它与 Lasso 的非零系数个数有什么区别？
  4. 从次梯度条件推导软阈值函数，并据此写出 Lasso 的 KKT 条件。
  5. 为什么只检查相邻两次系数变化不能保证坐标下降已经求得足够准确的解？
  6. 两个高度相关变量具有相近预测信息时，Ridge、Lasso 和 Elastic Net 可能分别怎样处理它们？
  7. 说明先对全样本标准化、再做交叉验证为什么会发生数据泄漏。哪些预处理步骤也有同样问题？
  8. 一标准误规则选择的模型为什么可能比最小交叉验证误差对应的模型更简单？它使用的“标准误”有什么局限？
  9. 区分系数稀疏和设计矩阵稀疏。一个问题能否只具有其中一种稀疏性？各举一例。
  10. 为什么稀疏 \(\mathbf X\) 的交叉乘积 \(\mathbf X^\top\mathbf X\) 仍可能是稠密矩阵？这对算法选择有什么影响？

14.11 上机实验（Lab）

  1. Lab 1：Ridge 的谱收缩。 人为构造具有指定奇异值的设计矩阵，比较 OLS 与不同 \(\lambda\) 下 Ridge系数对响应微小扰动的敏感程度，并把结果与 SVD 收缩因子联系起来。

  2. Lab 2：坐标下降诊断。 修改 coefficient_tol 和 kkt_tol，比较迭代次数、目标函数和 KKT 违反量。构造一个“系数变化很小但 KKT 仍不合格”的停止设置。

  3. Lab 3：暖启动。 分别用全零初值和上一 \(\lambda\) 的解计算完整 Lasso 路径，比较总迭代次数与最终系数差异。

  4. Lab 4：选择 \(\alpha\)。 在同一五折划分上比较\(\alpha\in\{0,0.25,0.5,0.75,1\}\)，对每个 \(\alpha\) 选择 \(\lambda\)，再在独立测试数据上评价最终模型。说明为什么不能为每个 \(\alpha\) 使用不同分折。

  5. Lab 5：选择稳定性。 改变变量组内相关程度和噪声标准差，重复 Bootstrap 选择频率实验，比较预测误差稳定性与变量集合稳定性。

  6. 拓展 Lab：稀疏存储。 逐步增大用户数、商品数和非零记录数，比较密集与稀疏表示的对象大小及矩阵—向量乘法时间。说明在哪些密度下稀疏表示不再具有明显优势。

14.12 延伸阅读

Ridge 的经典讨论见 Hoerl and Kennard (1970)，Lasso 见 Tibshirani (1996)，Elastic Net 见Zou and Hastie (2005)。Friedman, Hastie, and Tibshirani (2010) 系统说明了广义线性模型正则化路径的坐标下降算法；近端算法可参见 Parikh and Boyd (2014)，ADMM 可参见 Boyd et al. (2011)。

参考文献

Boyd, Stephen, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. 2011. “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.” Foundations and Trends in Machine Learning 3 (1): 1–122. https://doi.org/10.1561/2200000016.
Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. 2010. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” Journal of Statistical Software 33 (1): 1–22. https://doi.org/10.18637/jss.v033.i01.
Hoerl, Arthur E., and Robert W. Kennard. 1970. “Ridge Regression: Biased Estimation for Nonorthogonal Problems.” Technometrics 12 (1): 55–67. https://doi.org/10.1080/00401706.1970.10488634.
Parikh, Neal, and Stephen Boyd. 2014. “Proximal Algorithms.” Foundations and Trends in Optimization 1 (3): 123–231. https://web.stanford.edu/~boyd/papers/prox_algs.html.
Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society: Series B (Methodological) 58 (1): 267–88.
Zou, Hui, and Trevor Hastie. 2005. “Regularization and Variable Selection via the Elastic Net.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67 (2): 301–20. https://doi.org/10.1111/j.1467-9868.2005.00503.x.