R Матрица Гессе

Ответ на второй вопрос: «в основном да», и он предлагает почти немедленный ответ на ваш вопрос:

funcD <- sapply(vars, function(v) D(func, v))
funcDD <- matrix(list(), 2,2)
for (i in 1:length(vars)) 
        funcDD[,i] <- sapply(vars, function(v) D(funcD[[i]], v))
funcDD
#---------
     [,1]       [,2]      
[1,] Expression Expression
[2,] Expression Expression
> funcDD[1,1]
[[1]]
-(sin(x + y) + cos(x - y))

Квалификация «в основном» заключается в том, что нужно использовать «список», а не «выражение» в качестве типа объекта, который содержит матрица. Выражения на самом деле не квалифицируются как атомарные объекты, и вы можете легко извлечь значение и использовать его в качестве вызова, что может быть даже удобнее, чем иметь его в виде выражения:

> is.expression(funcDD[1,1])
[1] FALSE
> funcDD[1,1][[1]]
-(sin(x + y) + cos(x - y))
> class(funcDD[1,1][[1]])
[1] "call"

Оказывается, нужна была одна и та же структура, поэтому каждый элемент матрицы вызывается с тем же конкретным вектором, что и среда оценки, и возвращает их все в виде матрицы.:

matrix(sapply(funcDD, eval, env=list(x=0, y=pi)), length(vars))
#---------
     [,1] [,2]
[1,]    1   -1
[2,]   -1    1

Вот функция, которая может возвращать гессиан выражения в нескольких различных форматах. Код находится внизу этого ответа, ему предшествуют примеры его использования.

Примеры использования

my_fn <- expression((x^2)*(y^2))

# Get the symbolic Hessian as a character matrix

get_hessian(my_fn, as_matrix = TRUE)
#>      [x]              [y]             
#> [x] "2 * (y^2)"       "2 * x * (2 * y)"
#> [y] "2 * x * (2 * y)" "(x^2) * 2"

# Get the symbolic Hessian as a nested list of expressions
get_hessian(my_fn, as_matrix = FALSE)
#> $x
#> $x$x
#> 2 * (y^2)
#> 
#> $x$y
#> 2 * x * (2 * y)
#> 
#> 
#> $y
#> $y$x
#> 2 * x * (2 * y)
#> 
#> $y$y
#> (x^2) * 2

# Get the numeric Hessian from evaluating at a particular point
get_hessian(my_fn, eval_at = list(x = 2, y = 2))
#>      [x] [y]
#> [x]    8   16
#> [y]   16    8

Код для функции

get_hessian <- function(f, as_matrix = FALSE, eval_at = NULL) {

  fn_inputs <- all.vars(f); names(fn_inputs) <- fn_inputs
  n_inputs <- length(fn_inputs)

  # Obtain the symbolic Hessian as a nested list

  result <- lapply(fn_inputs, function(x) lapply(fn_inputs, function(x) NULL))

  for (i in seq_len(n_inputs)) {

    first_deriv <- D(f, fn_inputs[i])

    for (j in seq_len(n_inputs)) {

      second_partial_deriv <- D(first_deriv, fn_inputs[j])

      result[[i]][[j]] <- second_partial_deriv

    }
  }

  # Convert the symbolic Hessian to a character matrix
  if (is.null(eval_at)) {
    if (as_matrix) {
      matrix_result <- matrix(as.character(diag(n_inputs)), nrow = n_inputs, ncol = n_inputs)

      for (i in seq_len(n_inputs)) {
        for (j in seq_len(n_inputs)) {
          matrix_result[i, j] <- gsub("expression", "", format(result[[i]][[j]]), fixed = TRUE)
        }
      }

      dimnames(matrix_result) <- list(fn_inputs, fn_inputs)

      return(matrix_result)

    } else {

      return(result)

    }
  }

  # Evaluate the Hessian at a set point if a named list is provided

  if (!is.null(eval_at)) {
    result_vals <- diag(n_inputs)

    for (i in seq_len(n_inputs)) {
      for (j in seq_len(n_inputs)) {

        result_vals[i, j] <- eval(result[[i]][[j]], envir = eval_at)

      }
    }

    dimnames(matrix_result) <- list(fn_inputs, fn_inputs)

    return(result_vals)
  }
}

Вы можете использовать функцию hessian() из пакета calculus.

library(calculus)

# Create an expression with the function of interest
  func <- expression(sin(x+y)+cos(x-y))
  vars <- c("x", "y")

# Get the symbolic hessian
  hessian(f = func, var = vars)

# Get the hessian evaluated at a specific point
  hessian(f = func, var = c('x' = 0, 'y' = 1))

Я думаю, было бы намного проще написать цикл, который вычисляет каждую производную и помещает ее значение непосредственно в матрицу. Таким образом,

hess<-matrix(nrow=N,ncol=N)  #for x1 thru xN
for(j in 1:N) {
    for(k in 1:N) {
        hess[i,j]<- Dfunc(func,vars[i,j])
        }
    }

Где вам нужно будет настроить переменные x1,x2,...xN в матрице vars