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

Question

My aim is to call the hessian() function from the numDeriv R package from a cpp file (using Rcpp).

A toy example:
I want to calculate a hessian of a one-dimensional function x^n at the point x=1 with parameter n=3.
R code:

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

Cpp code:

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];

}

This works fine and I indeed get "6" at the output.
However my true goal is to calculate hessians of K-dimensional functions, therefore K=/=1. I try the following:

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

And in 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];

}

But now I get "invalid pointer" errors. A am not sure how to provide the hessian function call with a cpp function that takes a set of parameters to evaluate the partial derivatives at...

Answer 1

Couple of quick notes:

• Try the implementation in R and then move it to C++.
  • Provides a reference point and makes sure that everything works as intended.
• Search paths and names matter
  • Explicitly load numDeriv package before compiling.
  • Respect capitalization X vs. x.
• Ensure output types are accurate
  • From ?numDeriv::hessian, the output type is an N x N Rcpp::NumericMatrix instead of Rcpp::List.

Implementing in R

Coding the example and running it in pure R would give:

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.

So, immediately, this means that the k_dimensional() function is missing a reduction down to a scalar (e.g. single value).

Environment Run Time Error with C++ variant

After compiling the original code, there is a runtime error or when the function was called the issue an issue appears. For example, we have:

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

This provides the error of:

Error in call_KD_hessian_in_C(2) :

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

As we are using an R function found in a package not loaded by default, we must explicitly load it via library() or require() before calling the function.

Therefore, the process to avoid an environment issue should be:

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

Cleaned Up C++ Implementation

From prior discussion, note that we've:

1. Changed the function used with the hessian to be a scalar or single value, e.g. double, instead of a vector of values, e.g. NumericVector.
2. Ensured that before the function call the numDeriv R package is loaded.
3. Changed the return type expected from the hessian() function from Rcpp::List to Rcpp::NumericMatrix.

This results in the following C++ code:

#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;
}

Testing the routine gives:

# 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