c++ - 调用 numDeriv :hessian() with multiple-parameter-objective-function in Rcpp

Question

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

一个玩具示例:
我想计算参数 n=3 的一维函数 x^n 在点 x=1 处的粗麻布矩阵。
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];

}

但现在我收到“无效指针”错误。我不确定如何为 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++ 实现

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

1. 将与 hessian 一起使用的函数更改为标量 或单值，例如double，而不是值 vector ，例如NumericVector.
2. 确保在函数调用之前加载了 numDeriv R 包。
3. 将 hessian() 函数的预期返回类型从 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