Rcpp 中的环境

Environment in Rcpp

我正在处理涉及使用似然函数 (lik) 及其梯度 (grad) 的优化问题。这两个函数都使用 Rcpp 库进行编码,并提供给 R 中的 optim(...,method="BFGS") 优化函数。

在计算似然的过程中,创建了很多变量,对以后梯度的计算也很有用。为了避免每次需要梯度评估时都重新计算这些变量,我想出了创建一个环境(env)的想法,lik 使用 "communicate" 和 grad.换句话说,每次评估 lik 时,它都会将所需的变量 env 存储到环境中

likelihoodscript<-
using Rcpp::Environment; 
Environment e(env);
// likelihood function is computed
//var1 to varK are created
//they can be MatrixXd, SparseMatrix, double, etc.
e.assign("namevar1",var1);
e.assign("namevar2",var2);
...
e.assign("namevarK",varK);
return wrap(likvalue);

然后,当需要 grad 评估时,我从 grad 函数中搜索环境中所需的变量:

gradientscript<-
using Rcpp::Environment; 
Environment e(env);
//var1 to varK are recovered from the environment
SEXP input_var1 = e.get("var1");
int var1 = as<int>(var1);
SEXP input_var2= e.get("var2");
MatrixXd var2 = as<MatrixXd>(var2);
...
SEXP input_varK = e.get("varK");
double varK = as<double>(varK);
// gradient function is computed
return wrap(gradvector);

理论上,如果我使用相同的初始参数集开始优化问题,我每次都应该得到相同的结果(likgrad 中没有随机操作或 optim()), 但有时结果会有所不同。可以这么说,如果我用完全相同的参数重复优化 20 次,我有 15 次得到了很好的结果,但是其中 5 次我有一些不同(或非常不同)。我在这些错误中搜索了 "pattern",但它们看起来完全是随机的。

据我所知,这可能是由于在环境的散列 table 中错误的键值识别,因此我尝试在 R 中使用 hash 包。一些改进(也许18 倍的好结果,而不是 15),但还不是问题的确切解决方案。

如果有人知道问题出在哪里,请回答。 非常感谢。

P.D。我尝试了同样的事情,但使用 list 而不是 Environment:同样的问题,甚至更糟。

没有您身边的具体示例,很难判断问题出在哪里。这是我想出的可能是您的代码。也许这对您有帮助,或者它可以帮助您指出与我的示例不同的地方。 我的代码有效,所有不同的 optim 使用 cpp 函数运行都给出相同的结果。

library(inline)

#Rosenbrock Banana function as in optim examples
#it stores two values in the provided environment
fun_str='using Rcpp::Environment; 
Environment e(env);
NumericVector xx(x);
double a=(xx[1] - xx[0] * xx[0]);
double b=(1.0 - xx[0]);
e.assign("a",a);
e.assign("b",b);
return(wrap(100.0 * pow(a,2.0) + pow(b,2.0)));'

#corresponding gradient
#it accesses the two variables that we have stored above
grad_str='using Rcpp::Environment; 
Environment e(env);
double a=e["a"];
double b=e["b"];
NumericVector xx(x);
NumericVector res(2);
res[0]= -400.0 * xx[0] * a - 2.0 * b;
res[1]=200.0*a;
return(wrap(res));'

#define rcpp functions
fun_cpp=cxxfunction(signature(x="numeric",env = "environment"),
                    fun_str,
                    plugin = "Rcpp")
grad_cpp=cxxfunction(signature(x="numeric",env = "environment"),
                     grad_str,
                     plugin = "Rcpp")

res_vec=list()

#loop through 100 runs of the same BFGS optimization and store results
for (i in 1:100) {
  container_env=new.env()
  res_vec[[length(res_vec)+1]]=optim(par=c(-1.2,1),fn=fun_cpp,gr=grad_cpp,method="BFGS",env=container_env)
}

#compare all results
sum_var=0
for (i in 2:length(res_vec)) sum_var=sum_var+all.equal(res_vec[[i]],res_vec[[i-1]])
sum_var
#99

如您所见,它们都给出了相同的结果。 此外,结果与 optim 示例部分中给出的仅 R 代码的结果一致。