第 6 章优化算法

6.1 学习目标与路线

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

把最小二乘、极大似然和经验风险最小化写成目标函数，并识别参数定义域、约束与可导性；
从局部二次近似推导一元和多元 Newton 步，说明阻尼、线搜索与参数尺度为什么影响可靠性；
区分 Newton、BFGS、L-BFGS-B、Nelder–Mead、梯度下降和小批量 SGD 所使用的信息与适用条件；
核对解析梯度，读取 optim() 的收敛码、函数评价次数和最终梯度，而不是只抄录参数估计；
用 epoch、学习率衰减、验证损失和 checkpoint（某个 epoch 结束时保存的参数状态）组织小批量训练，并避免用测试集选择停止点或分类阈值；
结合目标函数值、尺度化梯度、对初值与尺度的敏感性以及必要的外部基准，判断计算结果是否可信。

阅读前建议熟悉第 3 章的浮点数与尺度、第 4 章的线性方程求解与条件数，以及第 5 章的 Newton 求根和失败诊断；同时需要梯度、Hessian 矩阵、最小二乘、极大似然和基本 R 矩阵运算。

本章沿着“问题设定—确定性优化—随机优化—诊断比较—完整案例”的路线展开。第五章已经说明 Newton法怎样求方程的根；本章只保留必要的桥接，并把重点放在下降方向、曲率、步长、约束和优化诊断上。Nelder–Mead 与近端思想帮助读者理解导数不可用或目标不可导时的处理方式；客户流失案例则把目标函数、预处理、训练、验证和测试连接成一个完整流程。

教学时可优先讲授优化问题的基本要素、Newton 与 BFGS、梯度下降、小批量 SGD、方法比较和客户流失案例。Nelder–Mead 与近端方法可作为选读内容；手写代码用于说明关键更新和诊断，不要求背诵。

6.2 优化问题的基本要素

上一章讨论的是求根问题：给定一个函数 \(g(x)\)，寻找使 \(g(x)=0\) 的位置。这个问题在统计计算中很重要，因为很多估计方程、分位数方程和置信区间端点都可以写成求根形式。本章要处理的问题更进一步：我们不只是要让某个函数等于 0，而是要在许多可能的参数中，找到使目标函数最大或最小的那一个。

在统计建模中，我们经常说“估计一个参数”“拟合一个模型”或“训练一个预测器”。从计算的角度看，这些说法背后常常是同一件事：在一组候选参数中，寻找使某个目标函数最好的一点。这个过程就是优化。一般地，优化问题可以写成

\[ \widehat{\boldsymbol\theta} =\arg\min_{\boldsymbol\theta\in\Theta} Q(\boldsymbol\theta), \]

其中 \(\boldsymbol\theta\) 表示待估计的参数向量，\(\Theta\) 表示参数允许取值的范围，\(Q(\boldsymbol\theta)\) 是我们希望尽量小的目标函数。如果目标是最大化某个函数，例如最大化对数似然\(\ell(\boldsymbol\theta)\)，也可以等价地写成最小化负对数似然 \(-\ell(\boldsymbol\theta)\)。

若参数空间 \(\Theta\) 为凸集，并且对任意可行点 \(\boldsymbol\theta_1,\boldsymbol\theta_2\) 和\(0\leq t\leq1\) 都有

\[ Q\{t\boldsymbol\theta_1+(1-t)\boldsymbol\theta_2\} \leq tQ(\boldsymbol\theta_1)+(1-t)Q(\boldsymbol\theta_2), \]

则称 \(Q\) 为凸函数。凸目标的任意局部极小点也是全局极小点；严格凸时至多有一个极小点。凸性仍不保证有限极小点一定存在，也不替代定义域和边界检查。

为什么统计计算需要优化算法？一个最直接的原因是，很多统计方法本来就是通过“最优”定义的。最小二乘估计是让残差平方和最小，极大似然估计是让观测数据出现的可能性最大，贝叶斯方法中的后验众数是让后验密度最大，机器学习中的经验风险最小化则是在训练集上寻找损失函数较小的模型。也就是说，优化不是某个额外的计算技巧，而是许多统计方法的基本语言。

在简单问题中，我们有时可以通过代数推导得到闭式解。例如，线性回归的最小二乘估计可以写成矩阵公式。但是，一旦模型稍微复杂，闭式解往往就不存在，或者即使存在也不适合直接计算。逻辑回归的极大似然估计通常没有显式公式，混合模型和隐变量模型的似然函数可能形状复杂，带有惩罚项的模型会引入不可导或约束条件，高维数据下还要考虑计算速度和数值稳定性。此时，我们不能只问“答案是什么公式”，还要问“怎样一步一步把答案算出来”。

这个联系也解释了为什么本章先从牛顿法讲起。上一章的牛顿法用于求解 \(g(x)=0\)；如果令\(g(x)=f'(x)\)，它就成为优化 \(f(x)\) 的方法。若一元函数 \(f(x)\) 可导，那么求 \(f(x)\) 的极值常常可以转化为求

\[ f'(x)=0. \]

在多元参数问题中，类似地，我们会寻找梯度为零的点。然而，统计建模中的目标函数不一定光滑，也不一定只有一个极值点；参数还可能受到非负性、概率和为 1、预算约束或正则化惩罚的限制。更重要的是，实际数据分析还关心初始值如何选、步长如何定、算法何时停止、结果是否稳定、是否只是局部最优。这些问题共同构成了优化算法的实践部分。

本章关注统计计算中最常用的优化思想。我们先从牛顿法出发，理解如何利用导数和曲率信息加快迭代；然后讨论多元优化问题中的梯度、海塞矩阵和收敛判据；接着介绍不依赖导数的 Nelder–Mead 算法，以及在大样本和机器学习中十分重要的梯度下降和随机梯度下降。最后，我们会通过比较不同算法的表现，并结合真实建模案例，理解优化算法在统计模型计算实现中的作用。

6.3 优化结果的基本检查

选择优化算法之前，需要先明确以下几个问题：

目标是什么？ 明确写出目标函数、参数单位以及最小化或最大化方向。对数似然通常取负后交给最小化程序；样本和、样本均值两种尺度虽然有同一个最优点，却会改变梯度和容忍误差的量级。
哪些参数是允许的？ 检查正值、概率、区间、线性等式等约束。能通过对数或 logit 变换自然满足的约束，可以考虑重新参数化；简单上下界可交给有边界的算法。
可利用什么结构？ 判断目标是否光滑、是否凸、梯度和海塞矩阵是否可得，以及目标能否拆成逐样本损失。算法应与这些结构相匹配。
数值尺度是否合理？ 参数量纲差异过大时，等高线会很狭长，停止准则也会失去可比性。标准化解释变量、设置参数尺度或采用合适的参数化，往往比盲目增加迭代次数更有效。
如何确认结果？ 至少报告终点的目标值、梯度或一阶最优性残差、迭代或函数评价次数、收敛状态和停止原因；非凸问题还应使用多个有意义的初值。若有闭式解或成熟软件，应把它作为基准。

所谓“算法收敛”只表示某个软件停止条件被触发，并不自动等于找到了全局最优点，也不表示模型具有良好的预测能力或因果解释。后文每个例子都会把优化输出与问题本身的诊断区分开来。

6.4 牛顿优化算法

6.4.1 从二次近似到 Newton 步

先考虑一元函数的无约束优化问题。为了叙述方便，本节主要讨论最小化问题：

\[ \min_x f(x). \]

如果要最大化 \(f(x)\)，可以等价地最小化 \(-f(x)\)。假设 \(f(x)\) 在感兴趣的区域内二阶可导，并且最优点位于区间内部。若 \(x^\ast\) 是局部极小点，那么通常需要满足一阶条件

\[ f'(x^\ast)=0. \]

这里要注意，“\(f'(x)=0\)”只是内部极值点的必要条件，不是充分条件。它也可能对应局部极大点、拐点，或者在多个驻点中只找到其中一个。因此，牛顿优化算法的核心不是简单地把导数设为 0，而是从一个初始值出发，利用局部导数信息逐步逼近一个合适的驻点。

牛顿法的更新公式可以从二次近似看出来。在当前点 \(x_n\) 附近，用二阶泰勒展开近似 \(f(x)\)：

\[ f(x)\approx f(x_n)+f'(x_n)(x-x_n)+\frac{1}{2}f''(x_n)(x-x_n)^2. \]

这个二次函数的极值点满足

\[ f'(x_n)+f''(x_n)(x-x_n)=0. \]

解出 \(x\)，就得到牛顿优化算法的迭代公式：

\[ x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}. \]

这个公式有一个直观解释：\(f'(x_n)\) 告诉我们当前位置的斜率，\(f''(x_n)\) 告诉我们函数在当前位置附近弯曲得有多厉害。牛顿法不是只沿着斜率方向走一步，而是利用曲率信息估计“如果局部形状近似为二次函数，极值点应该在哪里”。

这正是第 5.4 节对方程 \(g(x)=f'(x)=0\) 使用 Newton 求根。两章的关键差别是：这里还要确认驻点对应极小而非极大或拐点，并检查每一步是否下降。\(|f'(x_n)|\) 是一阶最优性监控量；参数或目标值几乎不再变化只能提示停滞，不能与小梯度共用一个“收敛”状态。若 \(f''(x_n)\) 接近 0 或为负，完整 Newton 步还可能很大或不是下降方向。

6.4.2 例：Poisson 均值的极大似然估计

下面用一个简单的极大似然问题说明牛顿优化算法如何出现在统计计算中，例如把 \(Y_i\) 看作企业第\(i\) 天收到的订单数。设\(Y_1,\ldots,Y_n\overset{\mathrm{iid}}{\sim}\operatorname{Poisson}(\lambda)\)，其观测值记为\(y_1,\ldots,y_n\)。为了保证 \(\lambda>0\)，令\(\eta=\log\lambda\)，即 \(\lambda=\exp(\eta)\)。忽略与参数无关的常数，负对数似然可以写成

\[ Q(\eta)=n\exp(\eta)-\left(\sum_{i=1}^n y_i\right)\eta. \]

对 \(\eta\) 求导得

\[ Q'(\eta)=n\exp(\eta)-\sum_{i=1}^n y_i,\qquad Q''(\eta)=n\exp(\eta). \]

因此，牛顿法可以用来最小化 \(Q(\eta)\)，从而得到 \(\lambda\) 的估计。

为了避免完整 Newton 步走得过远，下面使用 Armijo 回溯。令 \(p=-Q'(\eta)/Q''(\eta)\)；从\(\alpha=1\) 开始逐次减半，直到

\[ Q(\eta+\alpha p)\leq Q(\eta)+c_1\alpha Q'(\eta)p, \qquad 0<c_1<1. \]

右侧允许的下降量与当前位置的方向导数成比例，这称为“充分下降”条件。

set.seed(123)
y <- rpois(50, lambda = 3)

poisson_n <- length(y)
poisson_sum <- sum(y)
stopifnot(poisson_sum > 0)

poisson_objective <- function(eta) {
  poisson_n * exp(eta) - poisson_sum * eta
}
poisson_gradient <- function(eta) poisson_n * exp(eta) - poisson_sum
poisson_hessian <- function(eta) poisson_n * exp(eta)
poisson_scaled_gradient <- function(eta) {
  abs(poisson_gradient(eta)) /
    max(1, poisson_sum, poisson_n * exp(eta))
}

eta <- 0
poisson_trace <- data.frame(
    iter = 0,
    eta = eta,
    lambda = exp(eta),
    objective = poisson_objective(eta),
    scaled_gradient = poisson_scaled_gradient(eta),
    alpha = NA_real_,
    step = NA_real_
)
poisson_reason <- "达到最大迭代次数"

for (iter in seq_len(50)) {
    if (poisson_scaled_gradient(eta) <= 1e-10) {
      poisson_reason <- "尺度化梯度达到阈值"
      break
    }

    direction <- -poisson_gradient(eta) / poisson_hessian(eta)
    alpha <- 1
    current_value <- poisson_objective(eta)

    repeat {
      candidate <- eta + alpha * direction
      sufficient_decrease <- is.finite(poisson_objective(candidate)) &&
        poisson_objective(candidate) <= current_value +
          1e-4 * alpha * poisson_gradient(eta) * direction
      if (sufficient_decrease || alpha <= 2^-20) break
      alpha <- alpha / 2
    }

    if (!sufficient_decrease) {
      poisson_reason <- "线搜索未找到可接受步长"
      break
    }

    step <- alpha * direction
    eta <- eta + step
    poisson_trace <- rbind(
      poisson_trace,
      data.frame(
        iter = iter,
        eta = eta,
        lambda = exp(eta),
        objective = poisson_objective(eta),
        scaled_gradient = poisson_scaled_gradient(eta),
        alpha = alpha,
        step = step
      )
    )

    if (poisson_scaled_gradient(eta) <= 1e-10) {
      poisson_reason <- "尺度化梯度达到阈值"
      break
    }
}

poisson_fit <- list(
  eta = eta,
  lambda = exp(eta),
  objective = poisson_objective(eta),
  scaled_gradient = poisson_scaled_gradient(eta),
  converged = poisson_reason == "尺度化梯度达到阈值",
  reason = poisson_reason,
  trace = poisson_trace
)
knitr::kable(poisson_fit$trace, digits = 6, row.names = FALSE)

iter	eta	lambda	objective	scaled_gradient	alpha	step
0	0.000	1.000	50.00	0.677419	NA	NA
1	1.050	2.858	-19.87	0.078177	0.5	1.050000
2	1.135	3.111	-20.37	0.003399	1.0	0.084807
3	1.131	3.100	-20.37	0.000006	1.0	-0.003399
4	1.131	3.100	-20.37	0.000000	1.0	-0.000006


poisson_check <- data.frame(
  method = c("Newton", "closed form"),
  lambda = c(poisson_fit$lambda, mean(y))
)
knitr::kable(poisson_check, digits = 10, row.names = FALSE)

method	lambda
Newton	3.1
closed form	3.1

stopifnot(
  poisson_fit$converged,
  poisson_fit$scaled_gradient < 1e-8,
  abs(poisson_fit$lambda - mean(y)) < 1e-8
)

在这个例子中，\(\exp(\eta)\) 会逐渐接近样本均值。我们当然知道 Poisson 均值的极大似然估计有闭式解 \(\hat\lambda=\bar y\)，但这个例子有助于看清牛顿法的计算逻辑：每一步都用当前点的一阶导数和二阶导数修正参数。代码同时记录目标值和尺度化梯度，并用闭式解作外部核对；这比固定运行若干次后只查看参数更可靠。当模型更复杂、闭式解不存在时，同样的诊断原则仍然可以使用。

6.4.3 保护步与收敛诊断

牛顿法通常在最优点附近收敛很快，但它也有明显的限制。第一，它依赖初始值；初始值离目标点太远时，迭代可能走向不合适的驻点，甚至发散。第二，它需要二阶导数；在高维问题中，计算和求解海塞矩阵可能很昂贵。第三，如果目标函数不是凸函数，牛顿法找到的往往只是局部最优或某个驻点。因此，实际使用时通常要结合图形检查、多个初始值、步长控制和收敛诊断来判断结果是否可靠。

上面的实现没有盲目接受完整 Newton 步，而是先检查 Armijo 充分下降条件。如果候选点没有带来足够下降，就令步长系数 \(\alpha\) 依次取 \(1,1/2,1/4,\ldots\)。这种阻尼并不能把局部方法变成全局优化方法，却能避免很多由过大步长造成的失败。另一个重要区别是：梯度小反映一阶最优性，而步长小也可能来自曲率病态、边界或数值停滞；二者应当分别记录。若样本全为 0，上述 Poisson 模型的有限内部最优点并不存在，任何停止码也不能替代这一统计结构判断。

6.5 多元优化：Newton 法与拟 Newton 法

6.5.1 梯度、Hessian 矩阵与 Newton 方程

真实统计模型的参数通常不是一个数，而是一个向量。例如，线性回归中有截距和多个回归系数，逻辑回归中每个解释变量都有一个系数，混合模型和层次模型中还会同时包含均值、方差、权重等多类参数。因此，很多统计计算问题都可以写成下面的形式。为突出算法几何而不限定具体模型，本节把一般优化变量记为\(\mathbf{x}\)；当它表示模型参数时，可与前文的 \(\boldsymbol\theta\) 对应：

\[ \min_{\mathbf{x}\in\mathbb{R}^d} f(\mathbf{x}), \]

其中 \(\mathbf{x}=(x_1,\ldots,x_d)^\top\) 是 \(d\) 维优化变量，\(f(\mathbf{x})\) 是目标函数。与一元优化相比，多元优化的困难在于：我们不仅要判断“往左还是往右”，还要判断在高维空间中朝哪个方向移动。

多元优化中最重要的两个对象是梯度和海塞矩阵。梯度向量为

\[ \nabla f(\mathbf{x})= \left( \frac{\partial f(\mathbf{x})}{\partial x_1}, \frac{\partial f(\mathbf{x})}{\partial x_2}, \ldots, \frac{\partial f(\mathbf{x})}{\partial x_d} \right)^\top. \]

梯度可以理解为目标函数在当前位置上升最快的方向；如果要最小化函数，负梯度方向通常是首先想到的下降方向。海塞矩阵收集了所有二阶偏导数：

\[ \mathbf{H}(\mathbf{x}) =\left[ \frac{\partial^2 f(\mathbf{x})}{\partial x_i\partial x_j} \right]_{i,j=1}^d. \]

它描述的是目标函数在当前位置附近的弯曲程度。若 \(\mathbf{x}^\ast\) 是内部局部极小点，通常需要

\[ \nabla f(\mathbf{x}^\ast)=\mathbf{0}. \]

如果在该点附近海塞矩阵是正定的，即对任意非零向量 \(\mathbf{v}\) 都有\(\mathbf{v}^\top\mathbf{H}(\mathbf{x}^\ast)\mathbf{v}>0\)，那么目标函数在该点附近呈现“碗形”，这为局部极小提供了更强的证据。

多元牛顿法是上一节一元牛顿法的自然推广。在当前点 \(\mathbf{x}_n\) 附近，用二阶泰勒展开近似\(f(\mathbf{x})\)。令 \(\mathbf{p}=\mathbf{x}-\mathbf{x}_n\)，则有

\[ f(\mathbf{x}_n+\mathbf{p}) \approx f(\mathbf{x}_n) +\nabla f(\mathbf{x}_n)^\top\mathbf{p} +\frac{1}{2}\mathbf{p}^\top\mathbf{H}(\mathbf{x}_n)\mathbf{p}. \]

如果这个二次近似是一个凸的碗形函数，它的最小点满足

\[ \mathbf{H}(\mathbf{x}_n)\mathbf{p}_n=-\nabla f(\mathbf{x}_n). \]

于是我们得到更新公式

\[ \mathbf{x}_{n+1}=\mathbf{x}_n+\mathbf{p}_n. \]

这也是第 4 章线性代数方法在统计计算中的一个典型应用：实际编程时通常不显式计算\(\mathbf{H}^{-1}\)，而是求解线性方程组\(\mathbf{H}(\mathbf{x}_n)\mathbf{p}_n=-\nabla f(\mathbf{x}_n)\)。

带线搜索的多元 Newton 基本框架

选择初始值 \(\mathbf{x}_0\)，容忍误差 \(\epsilon\) 和最大迭代次数，设置 \(n=0\)。
计算梯度 \(\mathbf{g}_n=\nabla f(\mathbf{x}_n)\) 和海塞矩阵\(\mathbf{H}_n=\mathbf{H}(\mathbf{x}_n)\)。
求解线性方程组 \(\mathbf{H}_n\mathbf{p}_n=-\mathbf{g}_n\)。
检查 \(\mathbf{g}_n^\top\mathbf{p}_n<0\)。若不成立，修改 Hessian 矩阵或退回\(\mathbf{p}_n=-\mathbf{g}_n\)。
用 Armijo 回溯选择 \(0<\alpha_n\leq1\)，更新\(\mathbf{x}_{n+1}=\mathbf{x}_n+\alpha_n\mathbf{p}_n\)。
若 \(\lVert\mathbf{g}_{n+1}\rVert_2<\epsilon\)，报告满足一阶条件；若只有步长很小，则报告可能停滞；否则令 \(n=n+1\)，返回第 2 步。

常见的监控量包括：

\(\lVert\nabla f(\mathbf{x}_n)\rVert_2<\epsilon\)；
\(\lVert\mathbf{x}_n-\mathbf{x}_{n-1}\rVert_2<\epsilon\)；
\(|f(\mathbf{x}_n)-f(\mathbf{x}_{n-1})|<\epsilon\)。

后两项反映迭代是否还在移动，却不能单独证明一阶条件已经满足；小梯度也仍需结合曲率或凸性判断是极小点、极大点还是鞍点。因此实现时应区分“满足一阶条件”“停滞”和“确认极小点”。

6.5.2 例：二维凸目标函数

例6.1 用多元牛顿法最小化

\[ f(x_1,x_2)=x_1^2-x_1x_2+x_2^2+\exp(x_2). \]

这个函数的梯度和海塞矩阵分别为

\[ \nabla f(x_1,x_2)= \begin{pmatrix} 2x_1-x_2\\ -x_1+2x_2+\exp(x_2) \end{pmatrix}, \qquad \mathbf{H}(x_1,x_2)= \begin{pmatrix} 2 & -1\\ -1 & 2+\exp(x_2) \end{pmatrix}. \]

在这个例子中，海塞矩阵始终正定，因此目标函数只有一个局部极小点，也就是全局极小点。下面的 R代码从初始点 \((5,5)\) 出发，记录每次迭代的位置、目标函数值和梯度范数。

optim_f2 <- function(x) {
  x[1]^2 - x[1] * x[2] + x[2]^2 + exp(x[2])
}

optim_grad2 <- function(x) {
  c(2 * x[1] - x[2],
    -x[1] + 2 * x[2] + exp(x[2]))
}

optim_hess2 <- function(x) {
  matrix(c(2, -1,
           -1, 2 + exp(x[2])),
         nrow = 2, byrow = TRUE)
}

newton2d <- function(x0, tol = 1e-8, max_iter = 20) {
  if (length(x0) != 2 || any(!is.finite(x0))) {
    stop("x0 必须是长度为 2 的有限数值向量")
  }
  if (length(tol) != 1 || !is.finite(tol) || tol <= 0) {
    stop("tol 必须是有限正数")
  }
  if (length(max_iter) != 1 || !is.finite(max_iter) ||
      max_iter < 1 || max_iter != as.integer(max_iter)) {
    stop("max_iter 必须是正整数")
  }

  x <- x0
  if (!is.finite(optim_f2(x))) {
    stop("初始点的目标函数不是有限数")
  }
  path <- data.frame(
    iter = 0,
    x1 = x[1],
    x2 = x[2],
    value = optim_f2(x),
    grad_norm = sqrt(sum(optim_grad2(x)^2)),
    alpha = NA_real_
  )
  converged <- FALSE
  reason <- "达到最大迭代次数"

  for (iter in seq_len(max_iter)) {
    grad <- optim_grad2(x)
    if (sqrt(sum(grad^2)) <= tol) {
      converged <- TRUE
      reason <- "梯度范数达到阈值"
      break
    }

    hess <- optim_hess2(x)
    step <- solve(hess, -grad)
    if (!all(is.finite(step)) || sum(grad * step) >= 0) {
      reason <- "Newton 方程未给出有限下降方向"
      break
    }

    alpha <- 1
    current_value <- optim_f2(x)
    repeat {
      x_new <- x + alpha * step
      new_value <- optim_f2(x_new)
      sufficient_decrease <- is.finite(new_value) &&
        new_value <= current_value + 1e-4 * alpha * sum(grad * step)
      if (sufficient_decrease || alpha <= 2^-20) break
      alpha <- alpha / 2
    }
    if (!sufficient_decrease) {
      reason <- "线搜索未找到可接受步长"
      break
    }

    actual_step <- alpha * step
    x <- x_new
    grad_new <- optim_grad2(x)

    path <- rbind(
      path,
      data.frame(
        iter = iter,
        x1 = x[1],
        x2 = x[2],
        value = optim_f2(x),
        grad_norm = sqrt(sum(grad_new^2)),
        alpha = alpha
      )
    )

    if (sqrt(sum(grad_new^2)) <= tol) {
      converged <- TRUE
      reason <- "梯度范数达到阈值"
      break
    } else if (sqrt(sum(actual_step^2)) <=
               sqrt(.Machine$double.eps) * (1 + sqrt(sum(x^2)))) {
      reason <- "相对步长过小，可能发生数值停滞"
      break
    }
  }

  list(
    par = x,
    value = optim_f2(x),
    gradient_norm = sqrt(sum(optim_grad2(x)^2)),
    converged = converged,
    reason = reason,
    trace = path
  )
}

newton_fit <- newton2d(c(5, 5))
newton_path <- newton_fit$trace
round(newton_path, 4)
#>   iter      x1      x2    value grad_norm alpha
#> 1    0  5.0000  5.0000 173.4132  153.4946    NA
#> 2    1  1.9800  3.9600  64.2172   58.3961     1
#> 3    2  1.4388  2.8777  23.9840   22.0897     1
#> 4    3  0.8658  1.7316   7.8981    8.2467     1
#> 5    4  0.2890  0.5781   2.0332    2.6497     1
#> 6    5 -0.1146 -0.2291   0.8346    0.4515     1
#> 7    6 -0.2129 -0.4259   0.7892    0.0144     1
#> 8    7 -0.2163 -0.4326   0.7892    0.0000     1
#> 9    8 -0.2163 -0.4326   0.7892    0.0000     1
stopifnot(newton_fit$converged, newton_fit$gradient_norm < 1e-7)

图6.1: 二维牛顿法的迭代路径：红色点从初始值出发，逐步靠近目标函数的极小点。

从表格和图形可以看到，牛顿法前几步移动幅度较大，随后迅速稳定到极小点附近。它之所以快，是因为每一步都使用了梯度和海塞矩阵提供的局部形状信息。

不过，多元牛顿法也需要谨慎使用。若海塞矩阵接近奇异，方程组可能不稳定；若海塞矩阵不是正定的，牛顿方向甚至可能不是下降方向。实际软件中常会加入步长控制，例如使用

\[ \mathbf{x}_{n+1}=\mathbf{x}_n+\alpha_n\mathbf{p}_n,\qquad 0<\alpha_n\leq 1, \]

并通过线搜索选择合适的 \(\alpha_n\)。如果 \(\mathbf{g}_n^\top\mathbf{p}_n\geq 0\)，则当前 Newton方向不是下降方向，不能仅靠缩短步长修复；此时需要修改 Hessian 矩阵、退回负梯度方向，或改用拟Newton 方法。

6.5.3 梯度核对与 BFGS

在中等维度的光滑问题中，常用的默认选择不是显式计算 Hessian 矩阵，而是 BFGS 等拟 Newton 法。它利用相邻迭代的参数和梯度变化逐步更新曲率近似，通常比普通梯度下降少走很多“之”字形路径，又避免了每步形成和分解完整 Hessian 矩阵。有限内存版本 L-BFGS 只保存少量历史向量，适合参数更多的问题；R 的 optim() 用 L-BFGS-B 提供带简单上下界的版本。

解析梯度写错时，算法可能稳定地优化了一个与目标函数不一致的方向。下面用中心差分在若干参数点核对梯度。差分梯度不是更精确的“真值”，但作为独立实现，它能发现符号、常数因子和索引错误。

central_gradient <- function(fn, par,
                             rel_step = .Machine$double.eps^(1 / 3)) {
  step <- rel_step * pmax(1, abs(par))
  vapply(seq_along(par), function(j) {
    par_plus <- par
    par_minus <- par
    par_plus[j] <- par_plus[j] + step[j]
    par_minus[j] <- par_minus[j] - step[j]
    (fn(par_plus) - fn(par_minus)) / (2 * step[j])
  }, numeric(1))
}

check_points <- list(c(0.4, -0.6), c(2, 1), c(-1, -1.5))
gradient_check <- do.call(rbind, lapply(seq_along(check_points), function(i) {
  par <- check_points[[i]]
  analytic <- optim_grad2(par)
  numeric <- central_gradient(optim_f2, par)
  data.frame(
    point = i,
    coordinate = seq_along(par),
    analytic = analytic,
    finite_difference = numeric,
    absolute_difference = abs(analytic - numeric)
  )
}))
knitr::kable(gradient_check, digits = 8, row.names = FALSE)

point	coordinate	analytic	finite_difference
1	1	1.400	1.400
1	2	-1.051	-1.051
2	1	3.000	3.000
2	2	2.718	2.718
3	1	-0.500	-0.500
3	2	-1.777	-1.777

stopifnot(max(gradient_check$absolute_difference) < 1e-6)

中心差分要求正、负两个扰动点都可行；参数位于边界时，应改用单边差分或可行方向的方向导数。带有模拟噪声的目标以及尺度极端的参数还应检查多个步长，不能把一次吻合当作绝对证明。

通过核对后，再用 BFGS 求解同一个二维问题，并保留软件诊断。

bfgs_fit <- optim(
  par = c(5, 5),
  fn = optim_f2,
  gr = optim_grad2,
  method = "BFGS",
  control = list(reltol = 1e-12, maxit = 500)
)

bfgs_diagnostics <- data.frame(
  x1 = bfgs_fit$par[1],
  x2 = bfgs_fit$par[2],
  objective = bfgs_fit$value,
  gradient_norm = sqrt(sum(optim_grad2(bfgs_fit$par)^2)),
  function_evaluations = unname(bfgs_fit$counts["function"]),
  gradient_evaluations = unname(bfgs_fit$counts["gradient"]),
  convergence_code = bfgs_fit$convergence
)
knitr::kable(
  bfgs_diagnostics,
  digits = 8,
  col.names = c(
    "$x_1$", "$x_2$", "目标值", "梯度范数",
    "函数评价", "梯度评价", "收敛码"
  ),
  row.names = FALSE
)

\(x_1\)	\(x_2\)	目标值	梯度范数	函数评价	梯度评价	收敛码
-0.2163	-0.4326	0.7892	3.2e-07	22	12	0

bfgs_message <- if (is.null(bfgs_fit$message)) {
  "optim() 没有返回附加 message"
} else {
  bfgs_fit$message
}
cat(bfgs_message, "\n")
#> optim() 没有返回附加 message
stopifnot(
  bfgs_fit$convergence == 0,
  bfgs_diagnostics$gradient_norm < 1e-5
)

convergence == 0 表示 optim() 按其规则正常结束，不是“全局最优证明”。这里还检查了最终梯度，并可与前面的 Newton 结果互相核对。对于非凸目标，还应从多个有实际意义的初值重复运行，并比较终点目标值，而不是只保留最好看的一次结果。

6.5.4 边界约束与参数变换

若每个参数只有简单上下界，可以使用 L-BFGS-B。下面对同一目标增加 \(x_1\geq0\) 的约束。由于解落在边界，普通梯度不必等于 0；正确的一阶检查是：自由坐标的梯度接近 0，处在下界的坐标梯度不能指向可行域内的下降方向。

bounded_fit <- optim(
  par = c(0.5, 0),
  fn = optim_f2,
  gr = optim_grad2,
  method = "L-BFGS-B",
  lower = c(0, -Inf),
  upper = c(Inf, Inf),
  control = list(factr = 1e7, pgtol = 1e-10, maxit = 500)
)

bounded_gradient <- optim_grad2(bounded_fit$par)
projected_gradient <- bounded_gradient
at_lower <- abs(bounded_fit$par - c(0, -Inf)) < 1e-8
projected_gradient[at_lower & bounded_gradient > 0] <- 0

bounded_diagnostics <- data.frame(
  x1 = bounded_fit$par[1],
  x2 = bounded_fit$par[2],
  objective = bounded_fit$value,
  raw_gradient_norm = sqrt(sum(bounded_gradient^2)),
  projected_gradient_norm = sqrt(sum(projected_gradient^2)),
  convergence_code = bounded_fit$convergence
)
knitr::kable(bounded_diagnostics, digits = 8, row.names = FALSE)

x1	x2	objective	raw_gradient_norm	projected_gradient_norm	convergence_code
0	-0.3517	0.8272	0.3517	0	0

stopifnot(
  bounded_fit$convergence == 0,
  bounded_fit$par[1] >= -1e-10,
  bounded_diagnostics$projected_gradient_norm < 1e-5
)

正参数也常写成 \(\theta=\exp(\eta)\)，概率参数可写成\(p=\operatorname{logit}^{-1}(\eta)\)。重新参数化能让每个候选点天然可行，但会改变目标函数的曲率和先验或惩罚的含义；它不是纯粹的语法替换。一般线性或非线性约束还需要投影、拉格朗日乘子或专门的约束优化软件，本章不作系统展开。

6.6 无导数优化：Nelder–Mead 算法（选读）

前两节介绍的牛顿法很有效，但它依赖梯度和海塞矩阵。在实际统计计算中，我们有时只能计算目标函数的值，却很难写出导数。例如，目标函数可能来自一个模拟程序，可能包含排序、中位数、截断规则，也可能把数值积分、内层优化或外部黑箱程序包在里面。此时，仍然希望通过不断试探参数值来降低目标函数，这就需要不依赖导数的优化方法。

Nelder–Mead 算法就是最常用的无导数局部优化算法之一。R 中的 optim() 在没有特别指定方法时，默认使用的也是 Nelder–Mead 方法。需要注意的是，这里的“单纯形”不是线性规划中的单纯形法；它指的是由若干个点构成的几何形状。二维空间中的单纯形是三角形，三维空间中的单纯形是四面体；一般地，在 \(d\) 维空间中，一个单纯形由 \(d+1\) 个顶点构成。

图6.2: 二维单纯形示例：三个顶点构成一个三角形。

Nelder–Mead 算法维护一个单纯形，并在每一步比较这些顶点上的函数值。假设我们要最小化\(f(\mathbf{x})\)，当前单纯形的顶点按函数值从小到大排序为

\[ f(\mathbf{v}_1)\leq f(\mathbf{v}_2)\leq \cdots \leq f(\mathbf{v}_{d+1}). \]

其中 \(\mathbf{v}_1\) 是当前最好的点，\(\mathbf{v}_{d+1}\) 是当前最差的点。算法的基本想法很朴素：保留好的点，替换差的点，让单纯形逐渐向函数值更小的区域移动。

Nelder–Mead 算法

从 \(d+1\) 个初始顶点构造一个单纯形，并计算每个顶点的函数值。
按函数值排序，记最好的点为 \(\mathbf{v}_1\)，最差的点为 \(\mathbf{v}_{d+1}\)。
去掉最差点，计算其余 \(d\) 个点的中心

\[ \mathbf{c}=\frac{1}{d}\sum_{j=1}^{d}\mathbf{v}_j. \]

通过反射、扩张、收缩或压缩来产生新的单纯形。
若单纯形已经足够小，并且各顶点函数值差异已经足够小，则停止；否则返回第 2 步。

下面说明第 4 步中的几个基本动作。首先把最差点沿着中心 \(\mathbf{c}\) 的另一侧反射出去：

\[ \mathbf{v}_r=\mathbf{c}+\alpha(\mathbf{c}-\mathbf{v}_{d+1}),\qquad \alpha>0. \]

通常取 \(\alpha=1\)。记 \(f_r=f(\mathbf{v}_r)\)。标准决策还要把 \(f_r\) 同最好点 \(\mathbf{v}_1\)、次差点 \(\mathbf{v}_d\) 和最差点 \(\mathbf{v}_{d+1}\) 分别比较：若\(f(\mathbf{v}_1)\leq f_r<f(\mathbf{v}_d)\)，直接接受反射点；若 \(f_r<f(\mathbf{v}_1)\)，则沿该方向继续扩张：

\[ \mathbf{v}_e=\mathbf{c}+\gamma(\mathbf{v}_r-\mathbf{c}),\qquad \gamma>1. \]

常用取值是 \(\gamma=2\)，并在扩张点和反射点中保留目标值较小者。若反射点不优于次差点，则需要区分两种收缩。外收缩用于 \(f(\mathbf{v}_d)\leq f_r<f(\mathbf{v}_{d+1})\)：

\[ \mathbf{v}_{oc}=\mathbf{c}+\rho(\mathbf{v}_r-\mathbf{c}),\qquad 0<\rho<1. \]

只有 \(f(\mathbf{v}_{oc})\leq f_r\) 时才接受它。内收缩用于 \(f_r\geq f(\mathbf{v}_{d+1})\)：

\[ \mathbf{v}_{ic}=\mathbf{c}+\rho(\mathbf{v}_{d+1}-\mathbf{c}). \]

只有 \(f(\mathbf{v}_{ic})<f(\mathbf{v}_{d+1})\) 时才接受它。常用取值是 \(\rho=0.5\)。如果相应收缩仍然失败，说明当前单纯形可能太大或方向不合适，于是把除最好点以外的顶点都向最好点靠拢：

\[ \mathbf{v}_j \leftarrow \mathbf{v}_1+\sigma(\mathbf{v}_j-\mathbf{v}_1),\qquad j=2,\ldots,d+1,\quad 0<\sigma<1. \]

通常取 \(\sigma=0.5\)。这一步称为压缩，它会让搜索范围变小。

图 6.3 展示了 Nelder–Mead 算法优化函数\(f(x_1, x_2) = x_1^2 + x_2^2 + x_1\sin(x_2) + x_2\sin(x_1)\) 的过程。

图6.3: Nelder–Mead 单纯形迭代动图

图 6.3 展示了二维情形下单纯形如何逐步移动和缩小。和牛顿法不同，Nelder–Mead不需要知道梯度方向；它只是不断比较函数值，并根据“最好点”和“最差点”的相对位置调整单纯形。

实际使用时，停止条件通常不只看迭代次数，还会看两个方面：一是单纯形顶点之间的距离是否已经很小，二是各顶点对应的目标函数值是否已经非常接近。可以写成

\[ \max_{1\leq j\leq d+1}\lVert\mathbf{v}_j-\mathbf{v}_1\rVert_2<\epsilon_x, \qquad \max_{1\leq j\leq d+1}|f(\mathbf{v}_j)-f(\mathbf{v}_1)|<\epsilon_f. \]

这两个条件分别对应“参数位置变化不大”和“目标函数值变化不大”。本章生成图示的实现要求两者同时满足，避免在很平坦的方向上仅凭函数值接近而过早停止；其他软件可能采用不同的组合和尺度化规则。

下面用 optim() 演示 Nelder–Mead 算法。我们选择经典的 Rosenbrock 函数：

\[ f(x_1,x_2)=(1-x_1)^2+100(x_2-x_1^2)^2. \]

这个函数的全局最小点是 \((1,1)\)，但它的等高线呈狭长弯曲的谷地，因此常用来测试优化算法。

rosenbrock <- function(par) {
  x1 <- par[1]
  x2 <- par[2]
  (1 - x1)^2 + 100 * (x2 - x1^2)^2
}

nm_fit <- optim(
  par = c(-1.2, 1),
  fn = rosenbrock,
  method = "Nelder-Mead",
  control = list(maxit = 1000)
)

nm_fit$par
#> [1] 1.000 1.001
nm_fit$value
#> [1] 8.825e-08
nm_fit$convergence
#> [1] 0

convergence 等于 0 通常表示算法按软件设定的标准正常结束。这个例子也提醒我们，Nelder–Mead虽然使用方便，但并不保证找到全局最优。它对初始值和参数尺度比较敏感；如果不同参数的量纲差异很大，最好先做标准化，或在 optim() 中通过 control = list(parscale = ...) 调整参数尺度。对维度较高的问题，Nelder–Mead 往往会变慢，后面介绍的梯度下降和随机梯度下降更常用于大规模统计学习问题。

6.7 梯度方法

6.7.1 批量梯度下降

牛顿法利用二阶信息，Nelder–Mead 只利用函数值。梯度下降算法介于二者之间：它不需要海塞矩阵，但需要知道目标函数的一阶导数。对于许多统计模型，尤其是线性模型、逻辑回归和神经网络，目标函数的梯度相对容易计算，因此梯度下降成为统计计算和机器学习中非常重要的一类方法。

设 \(f(\mathbf{x})\) 是一个可微函数，我们希望最小化它。在当前位置 \(\mathbf{x}_n\)，梯度\(\nabla f(\mathbf{x}_n)\) 指向函数值上升最快的方向；因此，负梯度方向\(-\nabla f(\mathbf{x}_n)\) 是最自然的下降方向。梯度下降算法的基本更新为

\[\begin{equation} \mathbf{x}_{n+1}=\mathbf{x}_n-\alpha_n\nabla f(\mathbf{x}_n), \tag{6.1} \end{equation}\]

其中 \(\alpha_n>0\) 称为学习率或步长。学习率太小，算法会走得很慢；学习率太大，目标函数值可能来回震荡，甚至发散。一个理想化做法是在当前下降方向上寻找精确最优步长，即

\[\begin{equation} \alpha_n=\arg\min_{\alpha>0} f\left(\mathbf{x}_n-\alpha\nabla f(\mathbf{x}_n)\right). \tag{6.2} \end{equation}\]

这称为精确线搜索，但每一步都精确求解一个新的一维优化通常并不划算。更常用的是回溯线搜索：从\(\alpha=1\) 或某个试探值开始不断缩短步长，直到满足 Armijo 条件

\[ f(\mathbf{x}_n+\alpha\mathbf{p}_n) \leq f(\mathbf{x}_n)+c_1\alpha \nabla f(\mathbf{x}_n)^\top\mathbf{p}_n, \qquad 0<c_1<1. \]

这里 \(\mathbf{p}_n\) 必须是下降方向。实际应用中也会使用固定学习率，或者让学习率随着迭代逐渐变小。

梯度下降算法

选择初始值 \(\mathbf{x}_0\)、学习率规则和停止阈值，设置 \(n=0\)。
计算当前梯度 \(\nabla f(\mathbf{x}_n)\)。
按公式 (6.1) 更新参数。
若梯度范数足够小，报告满足一阶条件；若只有参数或目标变化很小，报告可能停滞；否则令\(n=n+1\)，返回第 2 步。

梯度下降的监控量可以参考牛顿法，例如\(\lVert\nabla f(\mathbf{x}_n)\rVert_2<\epsilon\)、\(\lVert\mathbf{x}_{n+1}-\mathbf{x}_n\rVert_2<\epsilon\)，或\(|f(\mathbf{x}_{n+1})-f(\mathbf{x}_n)|<\epsilon\)。后两项仍只能说明迭代变化小。还需要强调，只有在步长选择合适时，目标函数值才会稳定下降；公式本身并不自动保证每一步都变好。

图6.4: 一元函数中的梯度下降：从右侧初始点出发，沿负梯度方向逐步靠近极小点。

图6.5: 二元函数中的梯度下降：红色路径沿等高线的法线方向向中心移动。

在这个圆形等高线的例子中，负梯度方向几乎直接指向极小点，因此收敛很快。但实际统计模型中，参数之间常常存在相关关系，目标函数的等高线会变得狭长。此时负梯度方向可能在谷地两侧来回摆动，导致算法前进缓慢。这和第 3 章讨论的病态问题有相通之处：目标函数在不同方向上的曲率差异越大，优化越容易变慢。

图6.6: 当参数方向高度相关时，梯度下降路径可能出现明显的来回摆动。

例6.2 用梯度下降算法估计线性回归模型。

设线性回归模型为

\[ \mathbf{y}=\mathbf{X}\boldsymbol\beta+\boldsymbol\varepsilon. \]

其中 \(\mathbf y\in\mathbb R^n\) 是响应向量，\(\mathbf X\in\mathbb R^{n\times p}\) 是设计矩阵，\(\boldsymbol\beta\in\mathbb R^p\) 是参数向量，\(\boldsymbol\varepsilon\in\mathbb R^n\) 是误差向量。

给定 \(n\) 个观测样本，最小二乘损失可以写成

\[ J(\boldsymbol\beta) =\frac{1}{2n} (\mathbf{X}\boldsymbol\beta-\mathbf{y})^\top (\mathbf{X}\boldsymbol\beta-\mathbf{y}). \]

对应的梯度为

\[ \nabla J(\boldsymbol\beta) =\frac{1}{n}\mathbf{X}^\top (\mathbf{X}\boldsymbol\beta-\mathbf{y}). \]

下面先定义线性回归的损失函数、梯度和一个简单的梯度下降函数。为了让“自动选择步长”的例子运行更快，代码中使用了二次损失函数在线搜索方向上的解析最优步长。

lm.cost <- function(X, y, beta) {
  n <- length(y)
  sum((X %*% beta - y)^2) / (2 * n)
}

lm.cost.grad <- function(X, y, beta) {
  n <- length(y)
  as.vector(t(X) %*% (X %*% beta - y) / n)
}

lm.step.size <- function(X, y, beta) {
  grad <- lm.cost.grad(X, y, beta)
  Xgrad <- X %*% grad
  denom <- sum(Xgrad^2)
  if (denom <= .Machine$double.eps) return(0)
  length(y) * sum(grad^2) / denom
}

gd.lm <- function(X, y, beta.init, alpha, grad_tol = 1e-6,
                  max.iter = 100) {
  if (!is.matrix(X) || !is.numeric(X) || !is.numeric(y) ||
      !is.numeric(beta.init) || nrow(X) < 1 || nrow(X) != length(y) ||
      ncol(X) != length(beta.init) || any(!is.finite(c(X, y, beta.init)))) {
    stop("X、y 与 beta.init 必须维度匹配且全部有限")
  }
  if (!identical(alpha, "auto") &&
      (!is.numeric(alpha) || length(alpha) != 1 ||
       !is.finite(alpha) || alpha <= 0)) {
    stop("alpha 必须是 'auto' 或有限正数")
  }
  if (!is.numeric(grad_tol) || length(grad_tol) != 1 ||
      !is.finite(grad_tol) || grad_tol <= 0 ||
      !is.numeric(max.iter) || length(max.iter) != 1 ||
      !is.finite(max.iter) || max.iter < 1 ||
      max.iter != floor(max.iter)) {
    stop("grad_tol 必须为正，max.iter 必须是正整数")
  }

  use_line_search <- identical(alpha, "auto")
  beta <- beta.init
  betas <- list(beta)
  J <- list(lm.cost(X, y, beta))
  converged <- FALSE
  iter_done <- 0
  reason <- "达到最大迭代次数"

  for (iter in seq_len(max.iter)) {
    grad <- lm.cost.grad(X, y, beta)
    if (sqrt(sum(grad^2)) <= grad_tol) {
      converged <- TRUE
      reason <- "梯度范数达到阈值"
      break
    }

    step_alpha <- if (use_line_search) lm.step.size(X, y, beta) else alpha
    beta_new <- beta - step_alpha * grad

    if (!all(is.finite(beta_new))) {
      reason <- "迭代产生非有限参数"
      break
    }

    betas[[iter + 1]] <- beta_new
    J[[iter + 1]] <- lm.cost(X, y, beta_new)
    iter_done <- iter

    beta <- beta_new
  }

  final_grad_norm <- sqrt(sum(lm.cost.grad(X, y, beta)^2))
  if (final_grad_norm <= grad_tol) {
    converged <- TRUE
    reason <- "梯度范数达到阈值"
  }

  cat(reason, "\n")
  cat("共迭代", iter_done, "次\n")
  cat("系数估计为：", round(beta, 4), "\n")

  list(coef = betas, cost = J, niter = iter_done,
       gradient_norm = final_grad_norm,
       converged = converged, reason = reason)
}

set.seed(123)
beta0 <- 1
beta1 <- 3
sigma <- 1
n <- 10000
x <- rnorm(n, 0, 1)
y <- beta0 + x * beta1 + rnorm(n, mean = 0, sd = sigma)
X <- cbind(1, x)

gd.auto <- gd.lm(X, y, beta.init = c(-4, -5), alpha = "auto",
                 grad_tol = 1e-6, max.iter = 10000)
#> 梯度范数达到阈值 
#> 共迭代 2 次
#> 系数估计为： 0.9909 3.006

gd1 <- gd.lm(X, y, beta.init = c(-4, -5), alpha = 0.1,
             grad_tol = 1e-6, max.iter = 10000)
#> 梯度范数达到阈值 
#> 共迭代 154 次
#> 系数估计为： 0.9909 3.006

gd2 <- gd.lm(X, y, beta.init = c(-4, -5), alpha = 0.01,
             grad_tol = 1e-6, max.iter = 10000)
#> 梯度范数达到阈值 
#> 共迭代 1605 次
#> 系数估计为： 0.9909 3.006

ols_coef <- as.vector(qr.solve(X, y))
stopifnot(
  gd.auto$converged,
  gd1$converged,
  gd2$converged,
  max(abs(gd.auto$coef[[length(gd.auto$coef)]] - ols_coef)) < 1e-4
)

gd_summary <- data.frame(
  alpha = c("exact", "0.1", "0.01"),
  iterations = c(gd.auto$niter, gd1$niter, gd2$niter),
  objective = c(
    tail(unlist(gd.auto$cost), 1),
    tail(unlist(gd1$cost), 1),
    tail(unlist(gd2$cost), 1)
  ),
  gradient_norm = c(
    gd.auto$gradient_norm,
    gd1$gradient_norm,
    gd2$gradient_norm
  )
)
knitr::kable(
  gd_summary,
  digits = 7,
  caption = '不同步长规则下梯度下降的计算诊断',
  booktabs = TRUE,
  col.names = c("步长规则", "迭代次数", "目标函数", "梯度范数")
)

表6.1: 不同步长规则下梯度下降的计算诊断
步长规则	迭代次数	目标函数	梯度范数
exact	2	0.5015	5e-07
0.1	154	0.5015	9e-07
0.01	1605	0.5015	1e-06

从表 6.1 可以看到，步长会显著影响梯度下降的效率。自动线搜索通常能减少迭代次数；固定步长较小时，算法更稳，但移动较慢；固定步长过大时，算法则可能震荡或发散。图 6.7比较了两种固定步长下回归系数的收敛过程。

图6.7: 不同学习率下梯度下降算法的收敛过程。

6.7.2 小批量随机梯度下降

普通梯度下降每一步都要使用全部样本计算梯度。当样本量很大时，这一步可能非常昂贵。随机梯度下降（Stochastic Gradient Descent，SGD）的思想是：每次只用一个样本，或一小批样本，近似完整梯度，然后马上更新参数。以线性回归为例，完整梯度为

\[ \nabla J(\boldsymbol\beta)=\frac{1}{n}\sum_{i=1}^n \mathbf{x}_i(\mathbf{x}_i^\top\boldsymbol\beta-y_i). \]

SGD 在每次迭代中只抽取一个样本 \(i\)，用

\[ \nabla J_i(\boldsymbol\beta) =\mathbf{x}_i(\mathbf{x}_i^\top\boldsymbol\beta-y_i) \]

近似完整梯度。如果一次抽取 \(m\) 个样本，则称为小批量梯度下降，或 mini-batch gradient descent。

随机梯度下降算法

选择初始值 \(\boldsymbol\beta_0\)、学习率规则、批量大小 \(m\) 和最大 epoch 数。
在每个 epoch 开始时随机打乱样本，并把它们分成互不重叠的小批量。
依次用每个小批量计算近似梯度。
按照 \(\boldsymbol\beta\leftarrow \boldsymbol\beta-\alpha \widehat{\nabla J}(\boldsymbol\beta)\) 更新参数。
一个 epoch 结束后计算完整训练损失和必要的验证损失，保存 checkpoint。
重复第 2 到第 5 步，直到达到最大 epoch 数，或监控指标连续若干个 epoch 没有实质改善。

SGD 的单步计算成本低，特别适合大数据和在线学习。但它的路径通常带有随机波动，因为每一步看到的只是总体目标函数的一小部分。实际应用中常用的技巧包括：打乱数据顺序、使用 mini-batch、逐渐减小学习率，以及对变量做标准化。

一次小批量更新不等于一个 epoch。若样本量为 \(n\)、批量大小为 \(m\)，一个 epoch 大约包含\(\lceil n/m\rceil\) 次更新。下面的线性回归实现每个 epoch 无放回遍历全部样本，只在 epoch 末计算一次完整损失。这样既保留 SGD 的单步优势，也避免用偶然相近的两个随机批次判断“收敛”。

sgd.lm.cost <- function(X, y, beta) {
  n <- length(y)
  if (!is.matrix(X)) {
    X <- matrix(X, nrow = 1)
  }
  sum((X %*% beta - y)^2) / (2 * n)
}

sgd.lm.cost.grad <- function(X, y, beta) {
  n <- length(y)
  if (!is.matrix(X)) {
    X <- matrix(X, nrow = 1)
  }
  as.vector(t(X) %*% (X %*% beta - y) / n)
}

sgd.lm <- function(X, y, beta.init,
                   alpha = 0.1, batch_size = 256,
                   max_epochs = 150, patience = 15,
                   min_delta = 1e-8) {
  if (!is.matrix(X) || !is.numeric(X) || !is.numeric(y) ||
      !is.numeric(beta.init) || nrow(X) < 1 || nrow(X) != length(y) ||
      ncol(X) != length(beta.init) || any(!is.finite(c(X, y, beta.init)))) {
    stop("X、y 与 beta.init 必须维度匹配且全部有限")
  }
  controls <- c(alpha, batch_size, max_epochs, patience, min_delta)
  if (!is.numeric(controls) || length(controls) != 5 ||
      any(!is.finite(controls)) || alpha <= 0 || min_delta < 0 ||
      batch_size < 1 || batch_size != floor(batch_size) ||
      max_epochs < 1 || max_epochs != floor(max_epochs) ||
      patience < 1 || patience != floor(patience)) {
    stop("学习率与迭代控制参数不合法")
  }

  n <- length(y)
  beta <- beta.init
  checkpoints <- list(beta)
  train_loss <- sgd.lm.cost(X, y, beta)
  best_beta <- beta
  best_loss <- train_loss
  best_epoch <- 0
  epochs_without_improvement <- 0
  updates <- 0
  reason <- "达到最大 epoch 数"

  for (epoch in seq_len(max_epochs)) {
    shuffled_id <- sample.int(n)
    batch_id <- split(
      shuffled_id,
      ceiling(seq_along(shuffled_id) / batch_size)
    )
    epoch_rate <- alpha / sqrt(epoch)

    for (idx in batch_id) {
      grad <- sgd.lm.cost.grad(
        X[idx, , drop = FALSE],
        y[idx],
        beta
      )
      beta <- beta - epoch_rate * grad
      updates <- updates + 1
    }

    epoch_loss <- sgd.lm.cost(X, y, beta)
    checkpoints[[epoch + 1]] <- beta
    train_loss[epoch + 1] <- epoch_loss

    required_improvement <- min_delta * (1 + abs(best_loss))
    if (epoch_loss < best_loss - required_improvement) {
      best_loss <- epoch_loss
      best_beta <- beta
      best_epoch <- epoch
      epochs_without_improvement <- 0
    } else {
      epochs_without_improvement <- epochs_without_improvement + 1
    }

    if (epochs_without_improvement >= patience) {
      reason <- "完整训练损失在 patience 个 epoch 内没有实质改善"
      break
    }
  }

  gradient_norm <- sqrt(sum(sgd.lm.cost.grad(X, y, best_beta)^2))
  cat(reason, "\n")
  cat("完成", length(train_loss) - 1, "个 epoch，",
      updates, "次小批量更新\n")
  cat("最佳 checkpoint 系数：", round(best_beta, 4), "\n")

  list(
    par = best_beta,
    coef = checkpoints,
    cost = train_loss,
    niter = updates,
    epochs = length(train_loss) - 1,
    best_epoch = best_epoch,
    gradient_norm = gradient_norm,
    reason = reason
  )
}

现在在前面模拟的线性回归数据上运行小批量 SGD。参数轨迹会比普通梯度下降更有波动，但整体仍然朝着最小二乘解附近移动。

#> 完整训练损失在 patience 个 epoch 内没有实质改善 
#> 完成 29 个 epoch， 1160 次小批量更新
#> 最佳 checkpoint 系数： 0.9908 3.008

图6.8: 随机梯度下降算法的收敛过程。

SGD 的优势不在于每一步都比梯度下降更准确，而在于每一步更便宜。当样本量很大时，完整梯度下降可能要等很久才能更新一次参数；SGD 可以频繁更新，在相同计算预算下更早看到目标函数下降。它的代价是结果带有随机波动，因此通常需要配合学习率衰减、mini-batch 和多轮遍历数据来提高稳定性。这里的 patience 只表示完整训练损失暂时没有改善，不是精确的一阶最优证明；因此代码仍报告最佳checkpoint 的完整梯度范数，并用 QR 最小二乘解作核对。预测任务中更常监控独立验证集损失，案例部分将使用这一做法。

6.8 不可导目标与近端思想（选读）

前面介绍的 Newton 法、梯度下降和 SGD 都默认目标函数比较光滑，至少可以计算梯度。但在统计建模中，很多重要目标函数并不处处可导。例如，Lasso 的惩罚项包含 \(|\beta_j|\)，在 \(\beta_j=0\) 处不可导；约束优化也可以写成带有指标函数的优化问题，而指标函数通常不可导。对这类问题，与其强行求梯度，不如把目标函数拆成两个部分：

\[ \min_{\mathbf{x}} f(\mathbf{x})+g(\mathbf{x}), \]

其中 \(f\) 是光滑、容易求梯度的部分，\(g\) 是不可导但结构简单的部分。近端算法的核心思想是：对\(f\) 做普通梯度步，对 \(g\) 做一个“带二次惩罚的小优化”。给定步长 \(\alpha>0\)，函数 \(g\) 的近端算子定义为

\[\begin{equation} \operatorname{prox}_{\alpha g}(\mathbf{v}) = \arg\min_{\mathbf{x}} \left\{ g(\mathbf{x})+\frac{1}{2\alpha} \lVert \mathbf{x}-\mathbf{v}\rVert_2^2 \right\}. \tag{6.3} \end{equation}\]

这里的二次项要求新点不要离 \(\mathbf{v}\) 太远。如果 \(g\) 是某个集合的指标函数，近端算子就退化为向这个集合的投影；如果 \(g\) 是一范数惩罚，近端算子就是第 14 章 Lasso 小节会用到的软阈值函数。近端梯度法可以写成

\[\begin{equation} \mathbf{x}_{k+1} = \operatorname{prox}_{\alpha_k g} \left(\mathbf{x}_k-\alpha_k\nabla f(\mathbf{x}_k)\right). \tag{6.4} \end{equation}\]

公式 (6.4) 可以理解为“先沿着光滑部分下降一步，再用近端算子处理不可导部分”。在第 14 章的 Lasso 中，第一步处理最小二乘损失，第二步用软阈值函数把小系数压到 0。

标准凸收敛结论还需要条件：\(f\) 为凸函数且梯度是 \(L\)-Lipschitz 的，\(g\) 为适当的闭凸函数，步长通常取 \(0<\alpha\leq1/L\)。当 \(g\) 是集合的指标函数时，只有集合闭且凸，欧氏投影才对每个输入唯一。

最基本的例子是 \(g(x)=\lambda|x|\)。其近端算子有闭式解

\[ \operatorname{prox}_{\alpha\lambda|\cdot|}(v) =\operatorname{sign}(v)(|v|-\alpha\lambda)_+, \]

即软阈值：绝对值不超过阈值的数直接变成 0，其余数向 0 收缩。下面把闭式结果与一维数值优化核对。

soft_threshold <- function(v, threshold) {
  sign(v) * pmax(abs(v) - threshold, 0)
}

alpha <- 0.4
lambda <- 1.5
v_grid <- c(-2, -0.4, 0, 0.4, 2)
closed_form <- soft_threshold(v_grid, alpha * lambda)
numeric_check <- vapply(v_grid, function(v) {
  optimize(
    function(x) lambda * abs(x) + (x - v)^2 / (2 * alpha),
    interval = c(-4, 4),
    tol = 1e-12
  )$minimum
}, numeric(1))

prox_check <- data.frame(
  v = v_grid,
  closed_form = closed_form,
  numeric_optimization = numeric_check,
  absolute_difference = abs(closed_form - numeric_check)
)
knitr::kable(prox_check, digits = 8, row.names = FALSE)

v	closed_form	numeric_optimization
-2.0	-1.4	-1.4
-0.4	0.0	0.0
0.0	0.0	0.0
0.4	0.0	0.0
2.0	1.4	1.4

stopifnot(max(prox_check$absolute_difference) < 1e-6)

变量分裂和 ADMM 也能把复杂目标拆成容易处理的子问题，但它们需要同时诊断原始残差、对偶残差和惩罚参数。本节只建立近端直觉；第 14 章会在 Lasso 的具体更新中系统介绍这些内容。近端算法的进一步材料可见 Parikh and Boyd (2014)。

6.9 优化算法比较

到这里，我们已经看到几种不同类型的优化算法。它们没有绝对的优劣之分，关键是看问题本身：目标函数是否光滑？导数是否容易计算？样本量和参数维度有多大？我们是否更关心单步精确，还是更关心能否快速处理大规模数据？

算法	使用的信息	主要优点	主要限制	常见场景
牛顿法	梯度和海塞矩阵	在最优点附近通常收敛很快	需要二阶导数；海塞矩阵可能难算或不稳定	中低维、光滑、导数容易计算的问题
BFGS	目标函数和梯度	不显式计算 Hessian，通常是中等维光滑问题的可靠起点	需要保存稠密曲率近似；仍是局部方法	中等维极大似然和经验风险最小化
L-BFGS-B	目标函数、梯度和上下界	内存需求较低，可直接处理盒约束	不处理一般约束；边界处不能只看原始梯度	较高维光滑问题、参数有简单范围
Nelder–Mead	只使用函数值	不需要导数，适合黑箱目标函数	高维问题较慢；不保证全局最优	小维度、导数不可用的局部优化
梯度下降	梯度	实现简单，适合很多可微模型	对学习率敏感；病态问题中可能很慢	中大规模可微目标函数
小批量 SGD	小批量梯度	单次更新成本低，适合大数据和在线更新	路径有随机波动；需要调学习率并管理 checkpoint	大样本训练、流式数据
近端梯度	梯度和近端算子	利用“光滑损失+结构惩罚”	需要可快速计算的近端算子与合适步长	Lasso、稀疏与部分约束问题

这个表只是一个经验性的选择指南。实际建模时，我们还要考虑参数尺度、初始值、停止条件和数值稳定性。比较算法时不能只数外层迭代：一次 Newton 步需要 Hessian 和线性方程求解，一次 SGD 更新只看一个小批量。较公平的报告应同时给出最终目标值、最优性残差、函数和梯度评价次数、样本访问次数、运行时间与内存，并使用同一目标、起点和精度要求。

图6.9: 同一个线性回归损失函数下，牛顿法、梯度下降和随机梯度下降的迭代路径比较。

图 6.9 基于第 6.7.1 节模拟的线性回归数据，比较了三种可微优化算法在同一损失函数上的路径。这个例子是一个凸二次问题，对牛顿法非常有利；因此它不是所有优化问题的代表，而是用来帮助我们看清不同算法的计算风格。

牛顿法几乎直接到达最小点附近，因为线性回归的二次损失具有固定海塞矩阵。在这类问题中，二阶信息非常有价值。
梯度下降每一步只使用一阶信息，路径更长，但每一步的计算和实现都比较简单。它的表现很依赖步长和变量尺度。
小批量 SGD 的线只连接每个 epoch 的 checkpoint，星号单独标出训练损失最小的 checkpoint；两个相邻点之间已经包含多次批量更新。它牺牲单步精度，换取更低的单步计算成本。图上的点数因此不能直接当作计算成本。
Nelder–Mead 没有放在这张图里直接比较，因为它不利用可微结构。在小维度黑箱问题中它很有用；但对这种有明确梯度和海塞矩阵的二次损失，使用导数信息通常更高效。

在实际统计建模中，一个朴素但实用的选择顺序是：中等维光滑问题可先尝试带解析梯度的 BFGS，有简单上下界时使用 L-BFGS-B；二阶导数可靠且维度不高时可考虑带保护步的 Newton 法；目标可拆成大量逐样本损失时，再考虑小批量 SGD；只有梯度确实不可用的低维黑箱问题，才优先考虑 Nelder–Mead。无论选择哪种算法，都应检查目标函数值、最优性残差、停止信息和对初值与尺度的敏感性。

6.10 案例：客户流失预测——从优化到决策

前面几节分别介绍了牛顿法、Nelder–Mead、梯度下降和随机梯度下降。现在我们用一个客户流失预测问题把这些思想放到同一个建模流程中：给定客户的服务时长、费用、合同类型和支付方式等信息，预测客户是否已经流失。这里不仅要降低 Logistic 损失，还要用验证集选择训练 checkpoint 和分类阈值，最后把锁定的规则放到测试集评价。这样可以区分“优化目标下降”“样本外预测有效”和“决策规则合适”三个不同问题。

6.10.1 预测问题与 Logistic 目标函数

Logistic 回归常用于二分类问题。设二元随机变量 \(Y_i\in\{0,1\}\)，其观测值为 \(y_i\)；\(\mathbf{x}_i\) 是第 \(i\) 个样本的解释变量向量。模型假定

\[ P(Y_i=1\mid \mathbf{x}_i)=p_i,\qquad p_i=\frac{1}{1+\exp(-\eta_i)},\qquad \eta_i=\mathbf{x}_i^\top\boldsymbol\beta. \]

其中 \(p_i\) 是模型给出的流失概率，\(\eta_i\) 是线性预测值。在本节中，\(y_i=1\) 表示客户流失，\(y_i=0\) 表示客户未流失。把概率变成行动类别还需要阈值；阈值不是 Logistic 模型参数，应在模型训练之外依据验证数据和错误成本选择。

Logistic 回归通常通过极大似然估计参数。把极大化对数似然转化为最小化平均负对数似然，可写为

\[ J(\boldsymbol\beta)=\frac{1}{n}\sum_{i=1}^n \left\{\log(1+\exp(\eta_i))-y_i\eta_i\right\}. \]

令 \(\mathbf y=(y_1,\ldots,y_n)^\top\)、\(\mathbf p=(p_1,\ldots,p_n)^\top\)，并令设计矩阵\(\mathbf X\) 的第 \(i\) 行为 \(\mathbf x_i^\top\)，则梯度为

\[ \nabla J(\boldsymbol\beta) =\frac{1}{n}\mathbf{X}^\top(\mathbf{p}-\mathbf{y}). \]

这个形式与线性回归中的梯度非常相似，只是线性预测值先经过了 Logistic 函数。平均负对数似然始终是凸函数，因为其 Hessian 矩阵可写为\(\mathbf{X}^\top\mathbf{W}\mathbf{X}/n\)，其中 \(\mathbf{W}\) 为非负对角矩阵。完全分离不会破坏凸性，却可能使有限极小点不存在；设计矩阵的秩等条件还会影响严格凸性与解的唯一性。因此本例用 BFGS提供高精度确定性数值基准，再检查小批量 SGD 的实现。

6.10.2 稳定计算与 epoch 型 SGD

为了避免数值上溢，本节使用 plogis() 计算 Logistic 函数，并把单个观测的损失写成\(\log\{1+\exp[(1-2y)\eta]\}\)。这样可直接用稳定的 log1pexp() 计算，也避免在 \(y=1\) 且\(\eta\) 很大时出现“很大的数减很大的数”。

log1pexp <- function(z) {
  pmax(z, 0) + log1p(exp(-abs(z)))
}

logistic_loss <- function(X, y, beta) {
  eta <- drop(X %*% beta)
  mean(log1pexp((1 - 2 * y) * eta))
}

logistic_grad <- function(X, y, beta) {
  eta <- drop(X %*% beta)
  p <- plogis(eta)
  as.vector(crossprod(X, p - y) / length(y))
}

sgd.logisticReg <- function(X_train, y_train, X_valid, y_valid,
                            beta.init, alpha = 0.15,
                            batch_size = 128, max_epochs = 250,
                            patience = 20, min_delta = 1e-6) {
  if (!is.matrix(X_train) || !is.matrix(X_valid) ||
      !is.numeric(X_train) || !is.numeric(X_valid) ||
      !is.numeric(y_train) || !is.numeric(y_valid) ||
      !is.numeric(beta.init) ||
      nrow(X_train) < 1 || nrow(X_valid) < 1 ||
      ncol(X_train) != ncol(X_valid) ||
      nrow(X_train) != length(y_train) ||
      nrow(X_valid) != length(y_valid) ||
      length(beta.init) != ncol(X_train)) {
    stop("设计矩阵、响应变量与初始参数的维度不匹配")
  }
  if (any(!is.finite(X_train)) || any(!is.finite(X_valid)) ||
      any(!is.finite(beta.init)) || anyNA(y_train) || anyNA(y_valid) ||
      any(!y_train %in% c(0, 1)) || any(!y_valid %in% c(0, 1))) {
    stop("输入必须有限，响应变量必须取 0 或 1")
  }
  controls <- c(alpha, batch_size, max_epochs, patience, min_delta)
  if (!is.numeric(controls) || length(controls) != 5 ||
      any(!is.finite(controls)) || alpha <= 0 || min_delta < 0 ||
      batch_size < 1 || batch_size != floor(batch_size) ||
      max_epochs < 1 || max_epochs != floor(max_epochs) ||
      patience < 1 || patience != floor(patience)) {
    stop("学习率和迭代控制参数不合法")
  }

  n <- length(y_train)
  beta <- beta.init
  train_loss <- logistic_loss(X_train, y_train, beta)
  valid_loss <- logistic_loss(X_valid, y_valid, beta)
  trace <- data.frame(
    epoch = 0,
    train_loss = train_loss,
    valid_loss = valid_loss,
    train_gradient_norm = sqrt(sum(logistic_grad(X_train, y_train, beta)^2))
  )

  best_beta <- beta
  best_valid_loss <- valid_loss
  best_epoch <- 0
  patience_reference <- valid_loss
  epochs_without_improvement <- 0
  updates <- 0
  reason <- "达到最大 epoch 数"

  for (epoch in seq_len(max_epochs)) {
    shuffled_id <- sample.int(n)
    batches <- split(
      shuffled_id,
      ceiling(seq_along(shuffled_id) / batch_size)
    )
    epoch_rate <- alpha / sqrt(epoch)

    for (idx in batches) {
      batch_gradient <- logistic_grad(
        X_train[idx, , drop = FALSE],
        y_train[idx],
        beta
      )
      beta <- beta - epoch_rate * batch_gradient
      updates <- updates + 1
    }

    train_loss <- logistic_loss(X_train, y_train, beta)
    valid_loss <- logistic_loss(X_valid, y_valid, beta)
    trace <- rbind(
      trace,
      data.frame(
        epoch = epoch,
        train_loss = train_loss,
        valid_loss = valid_loss,
        train_gradient_norm = sqrt(sum(
          logistic_grad(X_train, y_train, beta)^2
        ))
      )
    )

    if (valid_loss < best_valid_loss) {
      best_beta <- beta
      best_valid_loss <- valid_loss
      best_epoch <- epoch
    }

    required_improvement <- min_delta * (1 + abs(patience_reference))
    if (valid_loss < patience_reference - required_improvement) {
      patience_reference <- valid_loss
      epochs_without_improvement <- 0
    } else {
      epochs_without_improvement <- epochs_without_improvement + 1
    }

    if (epochs_without_improvement >= patience) {
      reason <- "验证损失在 patience 个 epoch 内没有实质改善"
      break
    }
  }

  list(
    par = best_beta,
    trace = trace,
    best_epoch = best_epoch,
    best_valid_loss = best_valid_loss,
    train_loss = logistic_loss(X_train, y_train, best_beta),
    train_gradient_norm = sqrt(sum(
      logistic_grad(X_train, y_train, best_beta)^2
    )),
    epochs = max(trace$epoch),
    updates = updates,
    reason = reason
  )
}

函数每个 epoch 计算一次完整训练损失、验证损失和训练集梯度，并保存验证损失最小的参数。完整梯度在这里用于教学诊断，会额外遍历一次训练数据；更大规模的任务可以降低诊断频率。patience 表示允许连续多少个 epoch 没有实质改善；它既能容忍随机波动，也能避免无限训练。提前停止选中的是一个验证checkpoint，不保证训练集梯度已经接近 0，因此后面还会把完整梯度同 BFGS 基准比较。

6.10.3 训练、验证和测试数据

本案例使用 Hugging Face 上的 scikit-learn/churn-prediction 数据集（https://huggingface.co/datasets/scikit-learn/churn-prediction）。该数据记录了电信客户的账户信息和服务使用情况。原始数据中 TotalCharges 有少量空白值，本地文件data/telco_churn_hf.csv 保留了该变量可转为数值的 7032 条记录。

churn <- read.csv("data/telco_churn_hf.csv")
churn$churn01 <- as.integer(churn$Churn == "Yes")

churn_summary <- data.frame(
  指标 = c("样本量", "解释变量个数", "流失客户比例"),
  数值 = c(
    nrow(churn),
    8,
    round(mean(churn$churn01), 3)
  )
)
knitr::kable(churn_summary, caption = "客户流失数据的基本情况", row.names = FALSE)

表6.2: 客户流失数据的基本情况
指标	数值
样本量	7032.000
解释变量个数	8.000
流失客户比例	0.266

variable_table <- data.frame(
  变量 = c(
    "tenure", "MonthlyCharges", "TotalCharges", "SeniorCitizen",
    "Contract", "InternetService", "PaperlessBilling", "PaymentMethod", "Churn"
  ),
  含义 = c(
    "客户使用服务的月数",
    "每月费用",
    "累计费用",
    "是否为老年客户",
    "合同类型",
    "互联网服务类型",
    "是否使用无纸化账单",
    "支付方式",
    "是否流失"
  )
)
knitr::kable(variable_table, caption = "本节使用的变量", row.names = FALSE)

表6.3: 本节使用的变量
变量	含义
tenure	客户使用服务的月数
MonthlyCharges	每月费用
TotalCharges	累计费用
SeniorCitizen	是否为老年客户
Contract	合同类型
InternetService	互联网服务类型
PaperlessBilling	是否使用无纸化账单
PaymentMethod	支付方式
Churn	是否流失

这个数据同时包含数值变量和分类变量。连续变量需要标准化，以免尺度差异影响梯度方法；分类变量则用 model.matrix() 转换为 0/1 哑变量。我们按流失状态分层划分 60% 训练集、20% 验证集和 20%测试集。训练集用于更新参数，验证集用于选择 checkpoint 和阈值；测试集在规则锁定后才用于最终评价。

为了避免信息泄漏，均值和标准差只能从训练集估计，再把同一变换应用到验证集与测试集。哑变量仍保持 0/1，并没有被标准化；后面解释系数时必须保留这个差别。

set.seed(2024)

categorical_cols <- c("Contract", "InternetService", "PaperlessBilling", "PaymentMethod")
for (nm in categorical_cols) {
  churn[[nm]] <- factor(churn[[nm]], levels = sort(unique(churn[[nm]])))
}

stratified_three_way_split <- function(y, train_prop = 0.60,
                                       valid_prop = 0.20) {
  if (length(y) < 3 || anyNA(y)) {
    stop("y 必须至少包含 3 个非缺失观测")
  }
  if (!is.numeric(train_prop) || !is.numeric(valid_prop) ||
      length(train_prop) != 1 || length(valid_prop) != 1 ||
      !is.finite(train_prop) || !is.finite(valid_prop) ||
      train_prop <= 0 || valid_prop <= 0 ||
      train_prop + valid_prop >= 1) {
    stop("train_prop 与 valid_prop 必须为正，且两者之和小于 1")
  }
  train_id <- valid_id <- test_id <- integer(0)
  for (class_value in sort(unique(y))) {
    class_id <- sample(which(y == class_value))
    n_train <- floor(train_prop * length(class_id))
    n_valid <- floor(valid_prop * length(class_id))
    if (n_train < 1 || n_valid < 1 ||
        n_train + n_valid >= length(class_id)) {
      stop("每个类别在训练、验证和测试子集中都至少需要 1 个观测")
    }
    train_id <- c(train_id, class_id[seq_len(n_train)])
    valid_id <- c(valid_id, class_id[n_train + seq_len(n_valid)])
    test_start <- n_train + n_valid + 1
    test_id <- c(
      test_id,
      class_id[seq.int(test_start, length(class_id))]
    )
  }
  list(train = sort(train_id), valid = sort(valid_id), test = sort(test_id))
}

split_id <- stratified_three_way_split(churn$churn01)
train_data <- churn[split_id$train, ]
valid_data <- churn[split_id$valid, ]
test_data <- churn[split_id$test, ]

numeric_cols <- c("tenure", "MonthlyCharges", "TotalCharges")
train_center <- vapply(train_data[numeric_cols], mean, numeric(1))
train_scale <- vapply(train_data[numeric_cols], sd, numeric(1))
train_scale[train_scale == 0] <- 1

apply_training_scale <- function(data, center, scale) {
  transformed <- data
  transformed[names(center)] <- Map(
    function(x, mu, sig) (x - mu) / sig,
    data[names(center)], center, scale
  )
  transformed
}

train_data <- apply_training_scale(train_data, train_center, train_scale)
valid_data <- apply_training_scale(valid_data, train_center, train_scale)
test_data <- apply_training_scale(test_data, train_center, train_scale)

design_formula <- ~ tenure + MonthlyCharges + TotalCharges +
  SeniorCitizen + Contract + InternetService + PaperlessBilling + PaymentMethod

X_train <- model.matrix(design_formula, data = train_data)
y_train <- train_data$churn01
X_valid <- model.matrix(design_formula, data = valid_data)
y_valid <- valid_data$churn01
X_test <- model.matrix(design_formula, data = test_data)
y_test <- test_data$churn01

split_summary <- data.frame(
  subset = c("训练集", "验证集", "测试集"),
  n = c(length(y_train), length(y_valid), length(y_test)),
  churn_rate = c(mean(y_train), mean(y_valid), mean(y_test))
)
knitr::kable(
  split_summary,
  digits = 3,
  caption = "分层划分后的样本量与流失比例",
  row.names = FALSE
)

表6.4: 分层划分后的样本量与流失比例
subset	n	churn_rate
训练集	4218	0.266
验证集	1405	0.265
测试集	1409	0.266


stopifnot(
  length(unique(c(split_id$train, split_id$valid, split_id$test))) == nrow(churn),
  identical(colnames(X_train), colnames(X_valid)),
  identical(colnames(X_train), colnames(X_test)),
  all(is.finite(X_train)), all(is.finite(X_valid)), all(is.finite(X_test))
)

6.10.4 梯度核对、模型拟合与诊断

先在一个非零参数点用中心差分核对解析梯度，再分别运行 BFGS 和小批量 SGD。BFGS 使用全部训练样本寻找训练目标的极小点；SGD 则按验证损失保存最佳 checkpoint。两者回答的问题略有不同，但确定性解仍可帮助我们发现 SGD 代码中的严重错误。

set.seed(2025)
beta.init <- rep(0, ncol(X_train))

beta_check <- seq(-0.15, 0.15, length.out = ncol(X_train))
analytic_gradient <- logistic_grad(X_train, y_train, beta_check)
numeric_gradient <- central_gradient(
  function(beta) logistic_loss(X_train, y_train, beta),
  beta_check
)
logistic_gradient_check <- data.frame(
  max_absolute_difference = max(abs(analytic_gradient - numeric_gradient)),
  relative_l2_difference = sqrt(sum((analytic_gradient - numeric_gradient)^2)) /
    max(1, sqrt(sum(analytic_gradient^2)))
)
knitr::kable(logistic_gradient_check, digits = 9, row.names = FALSE)

max_absolute_difference	relative_l2_difference
0	0


bfgs_churn <- optim(
  par = beta.init,
  fn = function(beta) logistic_loss(X_train, y_train, beta),
  gr = function(beta) logistic_grad(X_train, y_train, beta),
  method = "BFGS",
  control = list(reltol = 1e-12, maxit = 1000)
)

res <- sgd.logisticReg(
  X_train, y_train, X_valid, y_valid,
  beta.init = beta.init,
  alpha = 0.8,
  batch_size = 128,
  max_epochs = 400,
  patience = 25,
  min_delta = 1e-6
)

beta_hat <- res$par
bfgs_gradient_norm <- sqrt(sum(
  logistic_grad(X_train, y_train, bfgs_churn$par)^2
))
optimization_diagnostics <- data.frame(
  method = c("BFGS：训练目标", "SGD：最佳验证 checkpoint"),
  train_loss = c(bfgs_churn$value, res$train_loss),
  valid_loss = c(
    logistic_loss(X_valid, y_valid, bfgs_churn$par),
    res$best_valid_loss
  ),
  train_gradient_norm = c(bfgs_gradient_norm, res$train_gradient_norm)
)
knitr::kable(
  optimization_diagnostics,
  digits = 6,
  caption = "BFGS 与小批量 SGD 的优化诊断",
  row.names = FALSE
)

表6.5: BFGS 与小批量 SGD 的优化诊断
method	train_loss	valid_loss	train_gradient_norm
BFGS：训练目标	0.4208	0.4314	0.000000
SGD：最佳验证 checkpoint	0.4209	0.4318	0.002158

work_diagnostics <- data.frame(
  method = c("BFGS", "SGD（含教学诊断）"),
  training_loss_evaluations = c(
    unname(bfgs_churn$counts["function"]),
    res$epochs + 2
  ),
  validation_loss_evaluations = c(1, res$epochs + 1),
  training_gradient_evaluations = c(
    unname(bfgs_churn$counts["gradient"]) + 1,
    res$epochs + 2
  ),
  epochs = c(NA, res$epochs),
  training_row_visits = c(
    length(y_train) * (sum(bfgs_churn$counts) + 1),
    length(y_train) * (3 * res$epochs + 4)
  )
)
knitr::kable(
  work_diagnostics,
  caption = "两种实现的计算工作量记录",
  col.names = c(
    "方法", "训练损失计算", "验证损失计算",
    "训练梯度计算", "Epoch", "训练样本行访问"
  ),
  row.names = FALSE
)

表6.6: 两种实现的计算工作量记录
方法	训练损失计算	验证损失计算	训练梯度计算	Epoch	训练样本行访问
BFGS	180	1	179	NA	1514262
SGD（含教学诊断）	213	212	213	211	2686866


stopifnot(
  logistic_gradient_check$max_absolute_difference < 1e-5,
  bfgs_churn$convergence == 0,
  bfgs_gradient_norm < 1e-5,
  all(is.finite(beta_hat)),
  res$best_valid_loss < res$trace$valid_loss[1]
)

“训练样本行访问”把一次完整训练损失、一次完整训练梯度和每个小批量更新所读取的训练行数相加；它没有计入验证集访问，也不等同于运行时间。表中 SGD 的教学实现为了每个 epoch 展示完整诊断而付出了明显监控成本，并在 211 个 epoch 内完成了 6963 次批量更新。计数包含返回最佳参数时的最终训练诊断；BFGS 计数包含 optim() 内部评价和表外的一次最终梯度检查，但不含拟合前的差分核对。这说明随机更新是否便宜还取决于检查频率。

图6.10: 小批量 SGD 的训练与验证损失；虚线表示验证损失最小的 checkpoint。

BFGS 给出高精度的全训练目标数值基准，而 SGD 保存的是验证损失最小的中间状态，所以 SGD 的训练梯度不一定像 BFGS 那样接近 0。这不是把一个差的停止码包装成“收敛”，而是明确区分数值优化终点与按样本外表现选择的 checkpoint。

只有 tenure、MonthlyCharges 和 TotalCharges 按训练集标准差进行了缩放；二元变量和哑变量仍是 0/1。因此不能把所有系数绝对值机械地排成统一的“变量重要性”。下表按设计矩阵顺序报告系数，并明确每个一单位变化的含义。系数表示控制表中其他变量后的预测关联，不是因果效应。

coef_table <- data.frame(
  variable = colnames(X_train),
  coefficient = as.numeric(beta_hat)
)
coef_table <- coef_table[coef_table$variable != "(Intercept)", ]
coef_table$unit <- ifelse(
  coef_table$variable %in% numeric_cols,
  "增加 1 个训练集标准差",
  "由 0 变为 1，或相对基准类别"
)
coef_labels <- c(
  tenure = "服务时长",
  MonthlyCharges = "每月费用",
  TotalCharges = "累计费用",
  SeniorCitizen = "老年客户",
  "ContractOne year" = "合同：一年",
  "ContractTwo year" = "合同：两年",
  "InternetServiceFiber optic" = "互联网服务：光纤",
  InternetServiceNo = "无互联网服务",
  PaperlessBillingYes = "无纸化账单",
  "PaymentMethodCredit card (automatic)" = "支付：信用卡自动扣款",
  "PaymentMethodElectronic check" = "支付：电子支票",
  "PaymentMethodMailed check" = "支付：邮寄支票"
)
coef_table$variable <- ifelse(
  coef_table$variable %in% names(coef_labels),
  unname(coef_labels[coef_table$variable]),
  coef_table$variable
)
coef_table$coefficient <- round(coef_table$coefficient, 3)

knitr::kable(
  coef_table,
  caption = "Logistic 回归系数及其一单位变化",
  col.names = c("变量", "系数", "一单位变化"),
  row.names = FALSE
)

表6.7: Logistic 回归系数及其一单位变化
变量	系数	一单位变化
服务时长	-1.208	增加 1 个训练集标准差
每月费用	0.082	增加 1 个训练集标准差
累计费用	0.403	增加 1 个训练集标准差
老年客户	0.318	由 0 变为 1，或相对基准类别
合同：一年	-0.623	由 0 变为 1，或相对基准类别
合同：两年	-1.500	由 0 变为 1，或相对基准类别
互联网服务：光纤	0.785	由 0 变为 1，或相对基准类别
无互联网服务	-0.721	由 0 变为 1，或相对基准类别
无纸化账单	0.403	由 0 变为 1，或相对基准类别
支付：信用卡自动扣款	-0.134	由 0 变为 1，或相对基准类别
支付：电子支票	0.398	由 0 变为 1，或相对基准类别
支付：邮寄支票	-0.069	由 0 变为 1，或相对基准类别

6.10.5 验证集选阈值与测试集评价

分类阈值应反映行动成本，而不是观察测试结果后随意挑选。为说明计算流程，假设漏判一个将流失客户的成本是误报一个未流失客户的 4 倍。这个比例只是教学设定，真实应用必须由业务后果、干预成本和容量共同确定。我们在验证集的候选阈值上最小化

\[ C(t)=1\times \operatorname{FP}(t)+4\times \operatorname{FN}(t), \]

然后锁定模型与阈值，再评价测试集。

prob_valid <- plogis(drop(X_valid %*% beta_hat))
threshold_grid <- seq(0.05, 0.95, by = 0.01)
threshold_diagnostics <- do.call(rbind, lapply(threshold_grid, function(t) {
  pred <- as.integer(prob_valid >= t)
  false_positive <- sum(pred == 1 & y_valid == 0)
  false_negative <- sum(pred == 0 & y_valid == 1)
  data.frame(
    threshold = t,
    false_positive = false_positive,
    false_negative = false_negative,
    cost = false_positive + 4 * false_negative
  )
}))
best_threshold_row <- threshold_diagnostics[
  which.min(threshold_diagnostics$cost),
]
threshold <- best_threshold_row$threshold

default_threshold_row <- threshold_diagnostics[
  which.min(abs(threshold_diagnostics$threshold - 0.5)),
]
threshold_comparison <- rbind(
  transform(best_threshold_row, rule = "验证集成本最小"),
  transform(default_threshold_row, rule = "固定阈值 0.5")
)
threshold_comparison$cost_per_100 <-
  100 * threshold_comparison$cost / length(y_valid)
knitr::kable(
  threshold_comparison[
    , c("rule", "threshold", "false_positive", "false_negative", "cost_per_100")
  ],
  digits = 2,
  caption = "验证集上的阈值选择",
  col.names = c("规则", "阈值", "误报", "漏判", "每百人加权成本"),
  row.names = FALSE
)

表6.8: 验证集上的阈值选择
规则	阈值	误报	漏判	每百人加权成本
验证集成本最小	0.2	369	51	40.78
固定阈值 0.5	0.5	109	181	59.29

有限验证集上的成本曲线可能很平，多个相邻阈值也可能并列。本例在给定网格中取第一个最小成本阈值；真实项目还应预先声明并列规则，并通过重采样考察阈值与成本的不确定性。

prob_test <- plogis(drop(X_test %*% beta_hat))
y_predicted <- as.integer(prob_test >= threshold)

confusion <- table(
  factor(y_predicted, levels = c(0, 1), labels = c("未流失", "流失")),
  factor(y_test, levels = c(0, 1), labels = c("未流失", "流失"))
)
names(dimnames(confusion)) <- c("预测值", "真实值")
knitr::kable(confusion, caption = "测试集混淆矩阵")

表6.9: 测试集混淆矩阵
	未流失	流失
未流失	639	43
流失	395	332

auc_rank <- function(prob, y) {
  if (!is.numeric(prob) || !is.numeric(y) || length(prob) == 0 ||
      length(prob) != length(y) || any(!is.finite(c(prob, y))) ||
      any(!y %in% c(0, 1))) {
    stop("prob 与 y 必须等长、有限，且 y 只能取 0 或 1")
  }
  r <- rank(prob, ties.method = "average")
  n1 <- sum(y == 1)
  n0 <- sum(y == 0)
  if (n1 == 0 || n0 == 0) return(NA_real_)
  (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

safe_ratio <- function(numerator, denominator) {
  if (denominator == 0) NA_real_ else numerator / denominator
}

true_positive <- unname(confusion["流失", "流失"])
true_negative <- unname(confusion["未流失", "未流失"])
false_positive <- unname(confusion["流失", "未流失"])
false_negative <- unname(confusion["未流失", "流失"])
accuracy <- mean(y_predicted == y_test)
sensitivity <- safe_ratio(true_positive, true_positive + false_negative)
specificity <- safe_ratio(true_negative, true_negative + false_positive)
precision <- safe_ratio(true_positive, true_positive + false_positive)
auc <- auc_rank(prob_test, y_test)
prob_clipped <- pmin(pmax(prob_test, 1e-15), 1 - 1e-15)
test_log_loss <- -mean(
  y_test * log(prob_clipped) + (1 - y_test) * log1p(-prob_clipped)
)
brier_score <- mean((prob_test - y_test)^2)
cost_per_100 <- 100 * (false_positive + 4 * false_negative) / length(y_test)

metric_table <- data.frame(
  metric = c(
    "验证集选择的阈值", "Log loss", "Brier 分数", "AUC", "准确率",
    "灵敏度", "特异度", "精确率", "每百人加权成本"
  ),
  value = c(
    threshold, test_log_loss, brier_score, auc, accuracy,
    sensitivity, specificity, precision, cost_per_100
  )
)
names(metric_table) <- c("指标", "数值")
knitr::kable(
  metric_table,
  digits = 3,
  caption = "锁定模型与阈值后的测试集表现",
  row.names = FALSE
)

表6.10: 锁定模型与阈值后的测试集表现
指标	数值
验证集选择的阈值	0.200
Log loss	0.419
Brier 分数	0.138
AUC	0.843
准确率	0.689
灵敏度	0.885
特异度	0.618
精确率	0.457
每百人加权成本	40.241

stopifnot(all(is.finite(metric_table$数值)))

Log loss 与训练目标一致，对把错误类别赋予极端概率的预测惩罚较重；Brier 分数是预测概率与 0/1结果的均方误差；AUC 衡量随机抽取的一对流失与未流失客户能否被正确排序，与某个固定阈值无关。这些指标关注的方面不同，任何一个都不能单独证明概率在所有人群中校准良好。

这个案例展示了三层检查。第一层是数值计算：梯度核对、BFGS 基准、训练损失和完整梯度；第二层是样本外预测：验证损失、测试 Log loss、Brier 分数和 AUC；第三层才是行动规则：在明确成本下选择阈值并观察误报与漏判。训练目标下降不能替代后两层证据。

模型系数和高风险预测只表示观察数据中的条件关联。它们不能说明某种合同或支付方式导致流失，也不能说明向高风险客户提供优惠一定产生更高的挽留收益。后一个问题还需要干预效果、客户价值与资源约束信息。优化器会忠实地优化输入给它的目标，却不会自动修复研究设计、概率校准或因果识别问题。

6.11 进一步阅读

关于线搜索、Newton 法、拟 Newton 法和约束优化，可参见 Nocedal and Wright (2006)；Bottou (2012) 讨论大规模学习中 SGD 的实践原则；Zeiler (2012) 从自适应学习率角度改进随机优化；Parikh and Boyd (2014) 系统介绍近端算子及其在 Lasso 等问题中的应用。阅读算法论文时，应特别区分理论迭代次数、函数或梯度评价次数、样本访问次数与实际运行时间。

6.12 本章小结

优化首先是问题设定，其次才是算法调用。目标函数、参数空间、约束、光滑性、可利用的逐样本结构和数值尺度共同决定算法选择。内部驻点需要小梯度，但小梯度不保证是最小值；边界最优点的原始梯度还可能不为 0。重新参数化与 L-BFGS-B 能处理常见的正值、概率和上下界约束，却会带来不同的曲率或边界诊断。

Newton 法利用 Hessian 矩阵，在合适的局部区域很快，但需要下降方向检查和线搜索等保护；BFGS 通过梯度变化近似曲率，是中等维光滑统计问题的实用起点；Nelder–Mead 适用于低维且导数确实不可用的黑箱目标。梯度下降对步长和变量尺度敏感，小批量 SGD 则用随机性换取低成本更新，需要用 epoch、shuffle、学习率计划、checkpoint 和验证监控来组织训练。不可导的结构化惩罚可以通过近端步骤处理。

可靠结果至少应报告目标值、最优性残差、函数和梯度评价或样本访问次数、停止原因、约束检查以及对初值和尺度的敏感性。解析梯度应先用独立差分核对，成熟实现或闭式解可作为基准。客户流失案例进一步说明，优化收敛、样本外预测和行动决策是三个层次：测试集不能参与 checkpoint 或阈值选择，预测关联也不能自动解释为干预效果。

6.13 思考题

为什么 \(\nabla Q(\boldsymbol\theta)=\mathbf{0}\) 只给出驻点候选，而不保证局部最小或全局最小？对\(Q(x)=(x-2)^2\) 且 \(0\leq x\leq1\)，最优点处的梯度为什么不为 0？
当 Hessian 矩阵不正定时，Newton 方向为什么可能使目标函数上升？回溯线搜索能修复过大的下降步，为什么不能修复一个本身不是下降方向的 Newton 步？
把收入变量从“元”改成“万元”不会改变回归模型所表达的预测关系，为什么会显著改变固定学习率梯度下降的轨迹和可用步长范围？
用中心差分核对梯度时，步长 \(h\) 过大和过小分别会产生什么误差？为什么一次参数点、一个 \(h\) 上的吻合仍不足以确认整个梯度函数正确？
optim() 返回 convergence = 0 后，还应检查哪些量？为什么多个初值到达同一点仍不是全局最优的严格证明？
比较 Newton、BFGS、梯度下降和 SGD 时，为什么只比较“迭代次数”不公平？一次迭代背后的计算工作可以用哪些指标描述？
随机优化中相邻两个小批量损失很接近，为什么可能是抽样巧合、学习率过小或停滞，而不是收敛？每次批量更新后计算完整训练损失又会怎样削弱 SGD 的优势？
为什么最佳验证 checkpoint 的训练梯度可以明显大于 BFGS 训练解的梯度？提前停止追求的是训练目标的一阶最优，还是样本外验证表现？
客户流失案例中，为什么不能把连续标准化变量与 0/1 哑变量的系数绝对值直接排成统一的重要性？
为什么分类阈值应依据事先声明的错误成本在验证集选择？“流失风险高”和“给予优惠会带来较高增量收益”之间还缺少哪些信息？

6.14 上机实验（Lab）

每份实验报告至少应包括目标函数与可行域、算法和初值、导数核对、关键控制参数、最终目标值与最优性残差、计算工作量、停止原因以及结果解释。随机算法还应报告种子、batch size、epoch、学习率规则和所选 checkpoint；失败案例同样要保留可复现证据。下列实验是题库，不要求一次完成全部内容；一学期教学可优先选择 Lab 1、Lab 2 和综合 Lab 5，其余按课时选做。

Lab 1：平滑目标的数值比较。 对 Rosenbrock 函数和本章二维凸函数，写出解析梯度；在多个参数点和多个差分步长下做方向导数核对。统一起点，比较手写梯度下降与 optim() 的 BFGS，报告目标值、梯度范数、函数和梯度评价次数、运行时间及停止信息；选做 Nelder–Mead 比较。改变起点后说明哪些结论稳定。
Lab 2：变量尺度、条件数与梯度下降。 模拟包含“年收入”和“失业率”等量纲差异明显变量的线性回归数据，用 QR 解作为基准。在原尺度和训练集标准化尺度运行固定步长与线搜索梯度下降，比较Hessian 条件数、可用学习率、迭代轨迹和梯度范数；把系数换回原单位并验证预测值一致。
Lab 3：边界、参数变换与最优性条件。 给定需求曲线\(D(p)=N\operatorname{logit}^{-1}(a-bp)\) 和单位成本 \(c\)，在 \(p\in[c,p_{\max}]\) 内最大化\((p-c)D(p)\)。比较 optimize() 与 L-BFGS-B，改变边界识别内部解和边界解，并检查投影梯度。再用 \(p=c+(p_{\max}-c)\operatorname{logit}^{-1}(\eta)\) 把无约束参数映射到价格区间内部，比较两种参数化的曲率与结果。说明有限 \(\eta\) 为什么取不到两个端点，以及真实边界最优时会出现什么数值现象；再解释观察性需求曲线为何不自动支持价格干预的因果解释。
Lab 4：小批量训练的计算预算。 在本章线性回归模拟中比较 batch size 为 1、32、256 和完整样本的更新。统一样本访问总量而不是更新次数，记录 epoch 损失、最佳 checkpoint、完整梯度范数和运行时间；再比较固定学习率与 \(\alpha_0/\sqrt{e}\) 衰减，解释随机性与收敛速度的折中。
综合 Lab 5：客户流失模型的优化与成本决策。 复用本章固定的数据划分，独立完成 Logistic 梯度核对、BFGS 基准与 epoch 型 SGD。比较至少两组学习率和 batch size 的训练/验证曲线，在验证集按两种误报—漏判成本选择模型 checkpoint 与阈值，锁定后报告测试 Log loss、Brier 分数、AUC、灵敏度、特异度和加权成本，并讨论预测关联、干预效果和客户价值的区别。
选做：黑箱与不可导目标。 第一部分在多个 parscale 和初值下用 Nelder–Mead 优化Rosenbrock 函数，并与有解析梯度的 BFGS 比较；第二部分实现软阈值与近端梯度，求解一个小型Lasso 问题并同第 14 章的成熟实现核对。分别说明无导数方法和近端方法利用或放弃了哪些目标结构。

参考文献

Bottou, Léon. 2012. “Stochastic Gradient Descent Tricks.” In Neural Networks: Tricks of the Trade, 421–36. Springer.

Nocedal, Jorge, and Stephen J. Wright. 2006. Numerical Optimization. 2nd ed. New York: Springer.

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.

Zeiler, Matthew D. 2012. “Adadelta: An Adaptive Learning Rate Method.” arXiv Preprint arXiv:1212.5701.

第 6 章 优化算法