在 Rcpp 中使用多参数 objective 函数调用 numDeriv:hessian()

Question

我的目标是从 cpp 文件（使用 Rcpp）调用 numDeriv R 包中的 hessian() 函数。

玩具示例：
我想计算参数 n=3.
在点 x=1 处的一维函数 x^n 的粗麻布矩阵 R代码：

H = call_1D_hessian_in_C(K=1)
print(H)

Cpp 代码：

double one_dimensional(double X, double N){
  return pow(X,N);
}

// [[Rcpp::export]]
double call_1D_hessian_in_C(double K) {

  Rcpp::Environment numDeriv("package:numDeriv");
  Rcpp::Function hessian = numDeriv["hessian"];
  double param = 3;

  Rcpp::List hessian_results =
  hessian(
    Rcpp::_["func"] = Rcpp::InternalFunction(one_dimensional), 
    Rcpp::_["x"] = 1.0,
    Rcpp::_["N"] = param
  );

  return hessian_results[0];

}

这很好用，我确实在输出中得到了“6”。
然而，我的真正目标是计算 K 维函数的粗麻布矩阵，因此 K=/=1。我尝试以下操作：

H = call_KD_hessian_in_C(K=2)
print(H)

在 Cpp 中：

NumericVector k_dimensional(NumericVector X, double N){
  return pow(X,N);
}

// [[Rcpp::export]]
double call_KD_hessian_in_C(double K) {

  Rcpp::Environment numDeriv("package:numDeriv");
  Rcpp::Function hessian = numDeriv["hessian"];
  double param = 3;

  Rcpp::NumericVector x = rep(1.0,K);

  Rcpp::List hessian_results = 
  hessian(
    Rcpp::_["func"] = Rcpp::InternalFunction(k_dimensional),
    Rcpp::_["x"] = x,
    Rcpp::_["N"] = param
  );

  return hessian_results[0];

}

但现在我收到 "invalid pointer" 错误。我不确定如何为 hessian 函数调用提供 cpp 函数，该函数采用一组参数来评估偏导数...

Answer 1

几个快速笔记：

尝试 R 中的实现，然后将其移至 C++。
- 提供参考点并确保一切按预期工作。
搜索路径和名称很重要
- 编译前显式加载 numDeriv 包。
- 尊重大写 X 与 x。
确保输出类型准确
- 从 ?numDeriv::hessian 开始，输出类型是 N x N Rcpp::NumericMatrix 而不是 Rcpp::List.

在 R 中实现

在 pure R 中对示例和运行进行编码将得到：

k = 2
k_dimensional = function(x, N) {
 x ^ N 
}

numDeriv::hessian(k_dimensional, x = rep(1, k), N = 3)

Error in hessian.default(k_dimensional, x = rep(1, k), N = 3) :

Richardson method for hessian assumes a scalar valued function.

所以，马上，这意味着 k_dimensional() 函数缺少对标量的缩减（例如单个值）。

环境运行 C++ 变体的时间错误

编译原始代码后，出现运行时错误或调用函数时出现问题。例如，我们有：

Rcpp::sourceCpp("path/to/call_KD_hessian_in_C.cpp")
call_KD_hessian_in_C(K = 2)

这提供了以下错误：

Error in call_KD_hessian_in_C(2) :

Cannot convert object to an environment: [type=character; target=ENVSXP].

由于我们使用的是在默认情况下未加载的包中找到的 R 函数，因此我们必须在调用之前通过 library() 或 require() 显式加载它函数。

因此，避免环境问题的过程应该是：

# Compile the routine
Rcpp::sourceCpp("path/to/call_KD_hessian_in_C.cpp")

# Load numDeriv to ensure it is on the search path
library("numDeriv")

# Call function
call_KD_hessian_in_C(2)

清理 C++ 实现

根据之前的讨论，请注意我们已经：

将与 hessian 一起使用的函数更改为标量或单值，例如double，而不是 值向量 ，例如NumericVector.
确保在函数调用之前加载 numDeriv R 包。
已将 hessian() 函数预期的 return 类型从 Rcpp::List 更改为 Rcpp::NumericMatrix。

这导致以下 C++ 代码：

#include <Rcpp.h>

double k_dimensional_cpp(Rcpp::NumericVector x, double N){
// ^^ Change return type from NumericVector to double

  // Speed up the process by writing the _C++_ loop
  // instead of relying on Rcpp sugar.
  double total = 0;
  for(int i = 0 ; i < x.size(); ++i) {
      total += std::pow(x[i], N);
  }

  // Return the accumulated total
  return total;
}

// [[Rcpp::export]]
Rcpp::NumericMatrix call_KD_hessian_in_C(double K) {

  // Ensure that numDeriv package is loaded prior to calling this function
  Rcpp::Environment numDeriv("package:numDeriv");
  Rcpp::Function hessian = numDeriv["hessian"];

  double param = 3;
  Rcpp::NumericVector x = Rcpp::rep(1.0, K);

  // Switched from Rcpp::List to Rcpp::NumericMatrix
  Rcpp::NumericMatrix hessian_results = 
  hessian(
    Rcpp::_["func"] = Rcpp::InternalFunction(k_dimensional_cpp),
    Rcpp::_["x"] = x,    // use lower case x to match function signature.
    Rcpp::_["N"] = param
  );

  // Return the calculated hessian
  return hessian_results;
}

测试例程给出：

# Ensure numDeriv is on search path
library("numDeriv")

# Call function
call_KD_hessian_in_C(K = 2)
#              [,1]         [,2]
# [1,] 6.000000e+00 3.162714e-12
# [2,] 3.162714e-12 6.000000e+00

在 Rcpp 中使用多参数 objective 函数调用 numDeriv:hessian()

Calling numDeriv:hessian() with multiple-parameter-objective-function in Rcpp

c++

statistics

r

linear-algebra

rcpp

在 R 中实现

环境运行 C++ 变体的时间错误

清理 C++ 实现