使用梯度函数 R 的数值 Hessian

Question

我需要使用我的梯度函数计算我的函数数值的Hessian（编程通过公式，不是数字）。 numDeriv 或 rootSolve 等包使用 numerical 梯度计算 hessian 不满足我的需要。我需要在内部执行的事情（我无法调用单独的方法）在 optim 包中，但是处理我的优化任务的唯一方法在 nlopt 包中很好地实现并将其最佳值传递给 optim 以获得 hessian 对我的程序来说太昂贵了。

所以我需要一些函数来计算 Hessian，而不是使用数字梯度（例如，参见这些公式 https://neos-guide.org/content/difference-approximations）。我无法创建这样的函数，因为我不明白如何选择参数 h（增量），我的函数对此非常敏感。我可以在 R 中找到这样的函数或以某种方式从 optim 包中检索它吗？或者至少有人可以解释如何为 h 选择错误最小化值，那么我将自己 post 这个函数？

Answer 1

如果我理解正确，你应该使用 numDeriv::jacobian()，它接受一个 vector-valued 函数并计算一个矩阵（每个元素相对于每个输入的导数），并将其应用于你的分析梯度函数。 jacobian() 确实使用了数值近似（准确地说是理查森外推法），但我没有看到任何其他方法可以从 black-box 梯度函数到 Hessian 矩阵？

您确实需要指定（或使用默认值）数值 'delta' 函数（默认为 1e-4）。另一方面，optim() 用于计算 Hessian 的内部代码也使用有限差分：参见 here and here ...

下面我定义了一个函数，它的梯度函数，还有它的Hessian；此代码显示 jacobian(grad(x)) 与 Hessian 相同（对于特定测试用例）。

library(numDeriv)
x1 <- c(1,1,1)

测试我没有搞砸梯度和 Hessian 函数：

all.equal(grad(f,x1),g(x1)) ## TRUE
all.equal(h(x1),hessian(f,x1)) ## TRUE

我的 Hessian 矩阵和梯度的 Jacobian 矩阵的数值等价性：

all.equal(h(x1),jacobian(g,x1)) ## TRUE

测试函数：

f <- function(x) {
  sin(x[1])*exp(x[2])*x[3]^2
}

g <- function(x) {
  c(cos(x[1])*exp(x[2])*x[3]^2,
    sin(x[1])*exp(x[2])*x[3]^2,
    sin(x[1])*exp(x[2])*2*x[3])
}  

h <- function(x) {
    m <- matrix(NA,3,3)
    m[lower.tri(m,diag=TRUE)] <-
        c(-sin(x[1])*exp(x[2])*x[3]^2,
          cos(x[1])*exp(x[2])*x[3]^2,
          cos(x[1])*exp(x[2])*2*x[3],
    # col 2
           sin(x[1])*exp(x[2])*x[3]^2,
           sin(x[1])*exp(x[3])*2*x[3],
    # col 3
           sin(x[1])*exp(x[2])*2)
    m <- Matrix::forceSymmetric(m,"L")
    m <- unname(as.matrix(m))
    return(m)
}

使用梯度函数 R 的数值 Hessian

Numeric Hessian using gradient function R

optimization

r

numeric

numerical-methods