Вызов numDeriv: hessian () с целевой функцией с несколькими параметрами в Rcpp - PullRequest
Новая
       47

Вызов numDeriv: hessian () с целевой функцией с несколькими параметрами в Rcpp

0 голосов
/ 20 марта 2019

Моя цель - вызвать функцию hessian() из пакета numDeriv R из файла cpp (используя Rcpp).

Пример игрушки:
Я хочу вычислить гессиан одномерной функции x ^ n в точке x = 1 с параметром n = 3.
R code: 

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

}

Но теперь я получаю ошибки "неверный указатель".Я не уверен, как обеспечить вызов функции гессиана с помощью функции cpp, которая принимает набор параметров для оценки частных производных при ...

1 Ответ

4 голосов
/ 20 марта 2019

Пара быстрых заметок:

  • Попробуйте реализацию в 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)

Ошибка в hessian.default (k_dimensional, x= rep (1, k), N = 3):

Метод Ричардсона для гессиана предполагает скалярную функцию.

Итак, это сразу означает, что функция k_dimensional()отсутствует редукция до скаляра (например, одно значение).

Ошибка времени выполнения среды с вариантом C ++

После компиляции исходного кода существует runtime error или при вызове функции возникает проблема.Например, у нас есть: 

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

Это обеспечивает ошибку:

Ошибка в call_KD_hessian_in_C (2):

Невозможно преобразовать объект в среду:[тип = символ;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. Изменена функция, используемая с гессианом, на скалярное или одиночное значение , например double, вместо вектора значения , например 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
...