Числового Гессиану с помощью градиента функции R

Question

мне нужно рассчитать гессенской моей функции численно с использованием моего градиент функции (программируется формулой, а не числовой). Пакеты, такие как numDeriv или rootSolve вычислить hessian с помощью численный градиент, который не удовлетворяет моим потребностям. Мне нужно выполнить внутренне (без отдельного метода, который я могу вызвать) в пакете optim, но единственный метод, который обрабатывает мою задачу оптимизации, хорошо выполненную в пакете nlopt и передавая ее оптимальное значение optim, чтобы получить hessian слишком дорого для моя программа.

Так что мне нужна функция, которая вычисляет гессиан, используя не числовой градиент (см., Например, эти формулы https://neos-guide.org/content/difference-approximations). Я не могу выполнить такую ​​функцию, поскольку не понимаю, как выбрать параметр h (increment), с которым моя функция очень чувствительна. Могу ли я найти такую ​​функцию в R или получить ее как-нибудь от optim? Или может кто-то хотя бы объяснить, как выбрать значение, минимизирующее ошибку, для h, поэтому я отправлю эту функцию самостоятельно?

Answer 1

Если вы правильно поняли, вы должны использовать numDeriv::jacobian(), который принимает векторнозначную функцию и вычисляет матрицу (производную от каждого элемента по отношению к каждому входу) и применяет ее к вашей аналитической градиентной функции. jacobian() действительно использует численное приближение (экстраполяция Ричардсона, если быть точным), но я не вижу другого способа, который вы могли бы получить от функции градиента черного ящика до гессиана?

Вам необходимо указать (или использовать значение по умолчанию) числовой функции «дельта» (по умолчанию 1e-4). С другой стороны, внутренний код, который optim() использует для вычисления гессиана также использует конечные различия: см here и here ...

Ниже я определить функцию, ее функцию градиента и его Гесс; этот код показывает, что jacobian(grad(x)) совпадает с гессианским (для конкретного тестового примера).

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

тест, что я не завинчивать градиент и Гессе функции:

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

Численный эквивалентность моей Гессе и якобиан градиента:

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) 
}