R- Hessian矩阵

Question

​

​

我需要创建一个函数的Hessian矩阵，如下所示：

func <- expression(sin(x+y)+cos(x-y))
vars <- c("x", "y")

我也需要二阶导数作为表达式，我需要评估它们很多次，所以我做了一个一阶导数的列表，和一个二阶导数列表。

funcD <- lapply(vars, function(v) D(func, v))
funcDD <- list(); for (i in 1:length(vars)) funcDD[[i]] <- lapply(vars, function(v) D(funcD[[i]], v))

到目前为止，它起作用了。

> funcDD
[[1]]
[[1]][[1]]
-(sin(x + y) + cos(x - y))

[[1]][[2]]
-(sin(x + y) - cos(x - y))


[[2]]
[[2]][[1]]
cos(x - y) - sin(x + y)

[[2]][[2]]
-(cos(x - y) + sin(x + y))

现在的问题是：，如何创建包含求值表达式值的矩阵？在外面试过了，没起作用。

> h <- outer(c(1:length(vars)), c(1:length(vars)), function(r, c) eval(funcDD[[r]][[c]], envir = list(x = 1, y = 2)))
Error in funcDD[[r]] : subscript out of bounds

其他问题：是否有一种更优雅的方法来存储二阶导数表达式？例如，是否可以将表达式存储在矩阵中而不是列表列表中？

第三个问题：可以得到表达式的变量向量吗？上面我使用了vars <- c("x"，"y")作为手动输入，它是必要的还是有"get_variables"-like方法？

Answer 1

第二个问题的答案是，“主要是肯定的”，它为你的问题提供了一个几乎立即的答案：

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

Answer 2

下面是一个函数，它可以以几种不同的格式返回表达式的Hessian。代码位于这个答案的底部，前面是它的使用示例。

示例用法

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

Answer 3

您可以使用来自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))