R中不同L-BFGS实现为何输出不同优化结果？

我在R中尝试用三种不同的L-BFGS实现优化同一个非凸函数：optim(method = "L-BFGS-B")、torch::optim_lbfgs()和lbfgs::lbfgs()。已统一设置最大迭代次数、梯度容差、函数容差、历史大小，使用相同起始点和目标函数，梯度计算统一采用numDeriv::grad()（Torch除外，用自动微分）。预期仅会因数值精度出现微小差异，但实际三者收敛到完全不同的最优解。复现代码如下：

library(torch)
library(lbfgs)
library(numDeriv)
library(ggplot2)

# Function to minimize
f <- function(x) sin(3 * x) * exp(-x^2) + 0.05 * x^2

# Starting points
starts <- c(-3, -2, 1.5)

# Common settings
max_iter <- 100
history_size <- 10
grad_tol <- 1e-8
f_tol <- 1e-9
factr_val <- f_tol / .Machine$double.eps  # For optim()

# Results table
finals <- data.frame(start = numeric(), method = character(),
                     x = numeric(), fx = numeric(), stringsAsFactors = FALSE)

for (x0 in starts) {
  ## --- optim (base R) ---
  res_optim <- optim(
    par = x0,
    fn = f,
    gr = function(x) numDeriv::grad(f, x),
    method = "L-BFGS-B",
    control = list(maxit = max_iter, lmm = history_size,
                   pgtol = grad_tol, factr = factr_val)
  )
  finals <- rbind(finals, data.frame(
    start = x0, method = "optim",
    x = res_optim$par, fx = res_optim$value
  ))

  ## --- torch ---
  x_torch <- torch_tensor(x0, requires_grad = TRUE)
  opt <- optim_lbfgs(
    list(x_torch), lr = 1, max_iter = max_iter,
    history_size = history_size, tolerance_grad = grad_tol,
    tolerance_change = f_tol, line_search_fn = "strong_wolfe"
  )
  closure <- function() {
    opt$zero_grad()
    loss <- torch_sin(3 * x_torch) * torch_exp(-x_torch^2) + 0.05 * x_torch^2
    loss$backward()
    loss
  }
  opt$step(closure)
  x_final <- as.numeric(x_torch)
  finals <- rbind(finals, data.frame(
    start = x0, method = "torch", x = x_final, fx = f(x_final)
  ))

  ## --- lbfgs (libLBFGS) ---
  res_lbfgs <- lbfgs(
    call_eval = f,
    call_grad = function(x) numDeriv::grad(f, x),
    vars = x0,
    m = history_size,
    epsilon = grad_tol,
    delta = f_tol,
    max_iterations = max_iter,
    linesearch_algorithm = "LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE"
  )
  x_final <- res_lbfgs$par
  finals <- rbind(finals, data.frame(
    start = x0, method = "lbfgs", x = x_final, fx = f(x_final)
  ))
}

# Print table
print(finals)

# Generate values
x_vals <- seq(-4, 4, length.out = 1000)
y_vals <- f(x_vals)
df <- data.frame(x = x_vals, y = y_vals)

# Plot
ggplot(df, aes(x = x, y = y)) +
  geom_line(color = "blue") +
  ggtitle("Function: sin(3x) * exp(-x^2) + 0.05 * x^2") +
  xlab("x") + ylab("f(x)") +
  theme_minimal()

核心差异原因 #

1. 线搜索实现的细节差异 #

尽管都指定了Strong Wolfe条件，但不同库的线搜索逻辑存在关键区别：

• optim()的L-BFGS-B采用More-Thuente线搜索的变种，其步长初始化、调整规则以及Wolfe条件的阈值判断和标准实现不同；
• torch::optim_lbfgs()的Strong Wolfe线搜索是PyTorch移植版本，步长更新的迭代次数限制、回溯系数等参数默认值与其他库不一致；
• lbfgs包基于libLBFGS库，其线搜索的最大迭代次数、c1/c2系数（Wolfe条件的参数）默认值与另外两个库不同，这些参数直接决定了每一步的步长选择，进而引导优化路径走向不同的局部极小值。

2. 梯度计算的精度与方式差异 #

• Torch使用自动微分，而optim()和lbfgs包使用numDeriv::grad()的数值微分。自动微分在浮点精度范围内是精确的，而数值微分存在截断误差，梯度的微小差异会在L-BFGS的迭代中被放大，导致搜索方向偏离；
• 即使同为数值微分，numDeriv::grad()的默认差分步长与libLBFGS内部可能隐含的梯度计算步长也可能不同，进一步加剧梯度差异。

3. 收敛停止条件的判断逻辑差异 #

虽然设置了相同的容差阈值，但不同库的判断逻辑并不等价：

• optim()的factr参数是通过f_tol / .Machine$double.eps转换的，实际判断的是目标函数的相对变化；而lbfgs包的delta参数是绝对函数容差，两者的停止条件触发时机不同；
• Torch的tolerance_change判断的是参数更新幅度，而非函数值或梯度的变化，与另外两个库的逻辑完全不同；
• 梯度容差的判断方式：有的库检查梯度的L2范数，有的检查最大分量的绝对值，这会导致优化器在不同迭代阶段停止，落入不同的局部极小。

4. L-BFGS核心更新规则的细节差异 #

不同库对L-BFGS的Hessian近似更新处理不同：

• 初始Hessian近似的选择：有的用单位矩阵，有的用最近一次更新的曲率信息缩放；
• 历史缓存的管理：当达到历史大小上限时，丢弃旧更新对的顺序和策略不同；
• 对非正定近似Hessian的修正策略：部分库会自动调整近似矩阵以保证正定性，而有的库则直接使用，这会影响搜索方向的计算。

验证与调试思路 #

• 统一线搜索参数：手动设置所有库的线搜索参数（如c1/c2系数、最大线搜索迭代次数、回溯系数），观察收敛结果是否更一致；
• 统一梯度计算方式：将Torch的自动微分替换为numDeriv::grad()，或让另外两个库使用自动微分（若可行），验证梯度差异的影响；
• 对齐收敛条件：修改参数让所有库使用相同的停止逻辑（如统一检查梯度L2范数小于阈值），对比收敛结果；
• 跟踪优化路径：记录每一步的参数值和目标函数值，对比三个优化器的搜索路径，定位差异出现的起始迭代步骤。

内容的提问来源于stack exchange，提问作者Alyde Bles