obtim_lbfgs Rcppnumerical 和 RcppEigen 中的函数与 Rcpparmadillo

Question

[编辑]

我正在寻找使用 Rcpparmadillo 的 Rcppnumerical 和 RcppEigen 中的 optim_lbfgs 函数的方法。我按照Rcpp Integration for Numerical Computing Libraries, but it was not working with the error, cannot declare variable 'obj' to be of abstract type 'scoreftn_mns'. But now, I fixed some codes to make it work. I define the beta as Eigen::VectorXd beta(p); as in the Rcppnumerical中的方式在模板中转换为arma::vec(?).

这是我正在尝试执行的代码。

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
#include <RcppNumerical.h>
using namespace Numer;
using namespace arma;
using namespace Rcpp;

// Eigen::Ref<Eigen::VectorXd>
// Eigen::Ref<const Eigen::VectorXd>

class socreftn_mns: public MFuncGrad
{
private:
  const arma::vec TIME;
  const arma::vec DELTA;
  const arma::mat COVARI;
  const arma::vec TARGETVEC;
public:
  socreftn_mns(const arma::vec Time, const arma::vec Delta, const arma::mat Covari,
               const arma::vec targetvec) : TIME(Time), DELTA(Delta), COVARI(Covari), TARGETVEC(targetvec) {}

  double f_grad(Constvec& beta, Refvec grad){

    arma::vec b_s = arma::vec(beta.data(),beta.size());

    int n = COVARI.n_rows;
    int p = COVARI.n_cols;

    arma::vec zero_vec_p = zeros(p);
    arma::mat zero_mat_np = zeros(n,p);

    arma::vec tempvec_p(p);
    arma::mat tempmat_np(n,p);

    arma::vec resid = log(TIME) + COVARI*b_s;

    arma::uvec index_resid = sort_index(resid);

    TIME(index_resid);
    DELTA(index_resid);
    COVARI.rows(index_resid);
    resid(index_resid);

    tempmat_np = zero_mat_np; arma::vec U_inf = zero_vec_p;
    for(int it=0; it<n; it++){
      tempmat_np = COVARI.row(it) - COVARI.each_row();
      U_inf += sum(tempmat_np.each_col()%conv_to<vec>::from((resid>=resid(it))),0).t()*DELTA(it);
    }
    U_inf = U_inf/n - TARGETVEC;

    double objvalue = conv_to<double>::from(sum(pow(U_inf,2)));

    double h = 1e-4;
    for(int itt=0; itt<p; itt++){
      tempvec_p = b_s;
      tempvec_p(itt) = tempvec_p(itt) + h;

      arma::vec resid_g = log(TIME) + COVARI*tempvec_p;

      arma::uvec index_resid_g = sort_index(resid_g);

      TIME(index_resid_g);
      DELTA(index_resid_g);
      COVARI.rows(index_resid_g);
      resid(index_resid_g);

      tempmat_np = zero_mat_np; arma::vec score_grad = zero_vec_p;
      for(int it=0; it<n; it++){
        tempmat_np = COVARI.row(it) - COVARI.each_row();
        score_grad += sum(tempmat_np.each_col()%conv_to<vec>::from((resid_g>=resid_g(it))),0).t()*DELTA(it);
      }
      score_grad = score_grad/n - TARGETVEC;

      double score_objvalue = conv_to<double>::from(sum(pow(score_grad,2)));

      grad(itt) = (score_objvalue-objvalue)/h;
    }

    // beta = Eigen::Ref(b_s);

    return objvalue;
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector aftsrr_bfgs(arma::vec Time, arma::vec Delta, arma::mat Covari, arma::vec targetvec){

  const arma::vec TIME = Time;
  const arma::vec DELTA = Delta;
  const arma::mat COVARI = Covari;
  const arma::vec TARGETVEC = targetvec;

  int p = COVARI.n_cols;
  // Score Function
  socreftn_mns obj(TIME, DELTA, COVARI, TARGETVEC);

  // Initial Guess
  Eigen::VectorXd beta(p);
  beta.setOnes();

  double fopt;
  int status = optim_lbfgs(obj, beta, fopt);
  if(status < 0)
    Rcpp::stop("fail to converge");

  return Rcpp::wrap(beta);
}

并且，这是有效的 R 代码。

library(Rcpp)
library(RcppArmadillo)
library(RcppEigen)
library(RcppNumerical)
library(survival)
library(aftgee)

sourceCpp("C:/Users/mattw/Documents/paper_wj/exercise/aftsrr_wj_cpp/aftsrr_wj/optim_bfgs.cpp")
U_beta_r_non = function(beta,Time,Delta,Covari) {

  n=length(Time)
  p=ncol(Covari)

  e_i_beta = as.vector(log(Time) + Covari %*% beta)

  order_resid = order(e_i_beta)

  Time = Time[order_resid]
  Covari = matrix(Covari[order_resid,],nrow = n)
  Delta = Delta[order_resid]
  e_i_beta = e_i_beta[order_resid]

  U_beta = list(NA)
  for (i in 1:n) {
    U_beta[[i]] = colSums(Delta[i]*t(Covari[i,]-t(Covari))*(e_i_beta>=e_i_beta[i]))
  }
  U_beta = Reduce('+',U_beta)/n
  U_beta = sum(U_beta^2)

  # grad = c(); h = 1e-4;
  # for (it in 1:p) {
  #   beta_g = beta;
  #   beta_g[it] = beta[it] + h
  #   
  #   e_i_beta = as.vector(log(Time) + Covari %*% beta_g)
  #   
  #   order_resid = order(e_i_beta)
  #   
  #   Time = Time[order_resid]
  #   Covari = matrix(Covari[order_resid,],nrow = n)
  #   Delta = Delta[order_resid]
  #   e_i_beta = e_i_beta[order_resid]
  #   
  #   grad_beta = list(NA);
  #   for (i in 1:n) {
  #     grad_beta[[i]] = colSums(Delta[i]*t(Covari[i,]-t(Covari))*(e_i_beta>=e_i_beta[i]))
  #   }
  #   grad_beta = Reduce('+',grad_beta)/n
  #   grad_beta = sum(grad_beta^2)
  #   
  #   grad[it] = (grad_beta-U_beta)/h
  # }

  # return(U_beta)
  return(sum(U_beta^2))
  # return(grad)
}

#------------------------DATA GENERATION----------------------

set.seed(1)

n=300
beta_0=1

gamma_0=0.5

Z1=matrix(rnorm(n,3,1),nrow=n)
Z2=matrix(rexp(n,5),nrow=n)

T_aft=as.vector(exp(-beta_0*Z1-gamma_0*Z2+rnorm(n,5,1)))
C_aft=as.vector(exp(-beta_0*Z1-gamma_0*Z2+rnorm(n,6,1)))
X_aft=C_aft*(T_aft>C_aft)+T_aft*(T_aft<=C_aft)
D_aft=0*(T_aft>C_aft)+1*(T_aft<=C_aft)
table(D_aft)

beta_aftsrr=-aftsrr(Surv(X_aft,D_aft)~Z1+Z2)$beta;beta_aftsrr

init_beta = rep(0,2)
cpp_bfgs_esti = aftsrr_bfgs(init_beta,X_aft,D_aft,cbind(Z1,Z2),rep(0,2));cpp_bfgs_esti
optim_lbfgsb = optim(init_beta,function(x){U_beta_r_non(x,X_aft,D_aft,cbind(Z1,Z2))},method = "L-BFGS-B")$par;optim_lbfgsb

现在，它正在工作，但是，它似乎只给出了 init_beta 而没有在 R 示例中进行任何计算。由于有人说上面的犰狳代码不能用，我正在研究Eigen库。但是，我不熟悉计算，很难转换它。所以，希望有一种方法可以通过一些修改来使用它。谢谢！

任何评论都会有所帮助。

Answer 1

我对编译您的代码时收到的其余错误消息感兴趣：

optim_num.cpp: In function ‘Rcpp::NumericVector aftsrr_bfgs(arma::vec, arma::vec, arma::mat, arma::vec)’:
optim_num.cpp:86:16: error: cannot declare variable ‘obj’ to be of abstract type ‘socreftn_mns’
   socreftn_mns obj(TIME, DELTA, COVARI, TARGETVEC);
                ^~~
optim_num.cpp:11:7: note:   because the following virtual functions are pure within ‘socreftn_mns’:
 class socreftn_mns: public MFuncGrad
       ^~~~~~~~~~~~
In file included from /usr/local/lib/R/site-library/RcppNumerical/include/integration/wrapper.h:13:0,
                 from /usr/local/lib/R/site-library/RcppNumerical/include/RcppNumerical.h:16,
                 from optim_num.cpp:6:
/usr/local/lib/R/site-library/RcppNumerical/include/integration/../Func.h:52:20: note:  virtual double Numer::MFuncGrad::f_grad(Numer::Constvec&, Numer::Refvec)
     virtual double f_grad(Constvec& x, Refvec grad) = 0;
                    ^~~~~~

从本质上讲，这意味着您的 class socreftn_mns 是从摘要 class MFuncGrad 派生的。这个 class 是抽象的，因为它不包含方法 f_grad(Constvec& x, Refvec grad) 的定义。你尝试通过定义方法f_grad(arma::vec& b_s, arma::vec grad)来定义这个，但是由于函数签名不同，虚函数没有被重载。因此你的 class 也是抽象的。

如果您使用相同的签名，应该可以解决问题。所需类型根据 Eigen 对象定义：

// Reference to a vector
typedef Eigen::Ref<Eigen::VectorXd>             Refvec;
typedef const Eigen::Ref<const Eigen::VectorXd> Constvec;

所以你必须在 Aramdillo 和 Eigen 结构之间来回转换。

obtim_lbfgs Rcppnumerical 和 RcppEigen 中的函数与 Rcpparmadillo

obtim_lbfgs function in Rcppnumerical and RcppEigen with Rcpparmadillo

optimization

rcpp

armadillo

r-package

numerical-analysis