Clash Royale CLAN TAG#URR8PPP

Calling Rcpp functions from optim

I'm trying to obtain MAP estimates for what basically amounts to a logistic regression model. I'm using the optim function which takes a log-posterior density along with its analytic gradient as arguments. I have both R versions and Rcpp versions of the density and gradient functions. I can successfully estimate the MAP with the R functions, but optim is entering asymtopia and failing to convergence to the optimum with the Rcpp functions.

I've verified that the R version of the density function and the Rcpp version of the density function return the same value:

ll_cpp = cpp_posterior_density(THETAi = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)

ll_R = lf_posterior_density(THETAi = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)


print(paste0(c("R: log posterior: ", ll_R)))
print(paste0(c("cpp: log posterior: ", ll_cpp)))


results in

"R: log posterior: " "15.8951804436067"
"cpp: log posterior: " "15.8951804436067"


I've also verified that the gradients are equal across the two versions.

d_cpp = grad(cpp_posterior_density, x = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)

d_R = grad(lf_posterior_density, x = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)

print(paste0(c("R: gradient of log posterior: ", paste(d_R, collapse = ", "))))
print(paste0(c("cpp: gradient of log posterior: ", paste(d_cpp, collapse = ", "))))


results in

[1] "R: gradient of log posterior: "
[2] "6.49720418347811, 4.67847452089852, 5.93682469664212, 1.47670777676947"

[1] "cpp: gradient of log posterior: "
"6.49720418347811, 4.67847452089852, 5.93682469664212, 1.47670777659075"


However, when I call optim with the Rcpp functions, I fail to get convergence :

#Using Rcpp
out_LF = optim(par = as.vector(THETA0_LF[i,]),
fn = cpp_posterior_density,
gr = cpp_grad_posterior_density,
Yi = as.vector(Y[i,]),
MUi = as.vector(MU_LF[i,]),
invS =invS,
TAU = TAU,
LAMBDA = LAMBDA,
J = J,
K = K,
method = "BFGS",
hessian = TRUE,
control = list(trace = 6)) #does not converge


results in

initial value 15.895180
final value -4748.586405


The final value must be strictly greater than zero indicating nonconvergence. However, with the R functions, I do get convergence:

#With R functions for density and gradient
out_LF2 = optim(par = as.vector(THETA0_LF[i,]),
fn = lf_posterior_density,
gr = lf_grad_posterior_density,
Yi = as.vector(Y[i,]),
MUi = as.vector(MU_LF[i,]),
invS =invS,
TAU = TAU,
LAMBDA = LAMBDA,
J = J,
K = K,
method = "BFGS",
hessian = TRUE,
control = list(trace = 6)) #converged


results in

initial value 15.895180
final value 11.980282


Any clue what is going on?

For reproducibility here is a link to a Dropbox folder with the required data (e.g., THETA0_LF, Y, MU_LF, etc.) and objective functions and gradients (both the R version and Rcpp version). Also included is an R file that reproduces the output above (see "debug-rcpp-for-credi.R").

Below is the Rcpp version of the objective function

#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]

double cpp_posterior_density(const arma::vec& THETAi, const arma::vec& Yi, const arma::vec& MUi, const arma::mat& invS, const arma::vec& TAU, const arma::mat& LAMBDA, const int J, const int K)

int j;
double lodd_j;
double b;
// PYi
arma::vec LT = LAMBDA*THETAi;
arma::vec PYi(J);
for (j = 0; j < J; j++)
lodd_j = LT(j) - TAU(j);
if(lodd_j<0)
b = 0;
else
b = lodd_j;

PYi(j) = exp(lodd_j-b)/(exp(-b) + exp(lodd_j-b));


double ll = 0.0;
for (j = 0; j < J; j++)
if (Yi(j)==1L)
ll += log(PYi(j));

if (Yi(j)==0L)
ll += log(1.0-PYi(j));



//Prior distriubtion
arma::vec dMUi = THETAi-MUi;
double twoprior = as_scalar(dMUi.t()*invS*dMUi);

// Return result
double dpost = -1.0*ll - 0.5*twoprior;
return dpost;



And below is the R version of the objective function:

lf_posterior_density<-function(THETAi, Yi, MUi, invS, TAU, LAMBDA,J,K, weight = NULL)

if (is.null(weight))weight = rep(1,J)

# Defined variables
# PYi - J (vector)
# ll - (scalar)
# dMUi - K (vector)
# prior - (scalar)

# Computations

PYi = as.vector(1/(1 + exp(TAU - LAMBDA%*%THETAi))) # J (vector)

# likelihood component
ll = as.numeric(0) #(scalar)
for (j in 1:J)
if (Yi[j] == 1L)ll = ll + weight[j]*log(PYi[j])
if (Yi[j] == 0L)ll = ll + weight[j]*log(1.0-PYi[j])


# prior distribution component
dMUi = (THETAi - MUi) # K (vector)
prior = as.numeric(-0.5*(dMUi%*%invS%*%dMUi)) #(scalar)

# Return
return(-ll - prior)




1 Answer
1

There is a difference in your objective functions:

ll_cpp = cpp_posterior_density(THETAi = 2*as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)

ll_R = lf_posterior_density(THETAi = 2*as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)

print(paste0(c("R: log posterior: ", ll_R)))
#> [1] "R: log posterior: " "22.495400131601"
print(paste0(c("cpp: log posterior: ", ll_cpp)))
#> [1] "cpp: log posterior: " "16.7463952181814"


I have not debugged your source code to find the error.

In such cases it is useful to add REPORT = 1 to the control list. For R this gives:

REPORT = 1

control

initial value 45.707620
iter 2 value 28.881100
iter 3 value 22.426070
iter 4 value 20.145499
iter 5 value 19.922129
iter 6 value 19.805083
iter 7 value 19.684769
iter 8 value 19.684366
iter 9 value 19.684345
iter 10 value 19.684343
iter 10 value 19.684343
final value 19.684343
converged


For Rcpp:

initial value 45.707620
iter 2 value 23.059207
iter 3 value -33.279972
iter 4 value -77.878965
iter 4 value -77.878965
iter 5 value -93.872445
iter 5 value -93.872445
iter 6 value -2830.594586
iter 6 value -2830.594586
iter 6 value -2830.594586
final value -2830.594586
converged


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.