3

黑森矩阵

我需要创建一个函数的 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”的方法？

4

4 回答 4

4

问题二的答案是“大部分是的”，它几乎可以立即回答您的问题：

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
于 2013-11-30T15:29:30.097 回答
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)
  }
}
于 2019-08-19T02:17:28.047 回答
1

您可以使用包中的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))
于 2020-01-08T19:48:31.260 回答
0

我认为编写一个循环来计算每个导数并将其值直接放入矩阵中会容易得多。因此，

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

于 2013-11-30T13:33:27.593 回答