带有 Rcpparmadillo 的 Rcppnumerical 和 RcppEigen 中的 obtim_lbfgs 函数

Posted 2023-02-16

【中文标题】带有 Rcpparmadillo 的 Rcppnumerical 和 RcppEigen 中的 obtim_lbfgs 函数

【英文标题】：obtim_lbfgs function in Rcppnumerical and RcppEigen with Rcpparmadillo

【发布时间】：2018-06-07 15:36:39

【问题描述】：

[编辑]

我正在寻找将 Rcppnumerical 和 RcppEigen 中的 optim_lbfgs 函数与 Rcpparmadillo 一起使用的方法。我按照Rcpp Integration for Numerical Computing Libraries 中的方式进行操作，但它无法解决错误cannot declare variable 'obj' to be of abstract type 'scoreftn_mns'。但是现在，我修复了一些代码以使其正常工作。我将beta 定义为Eigen::VectorXd beta(p);，就像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 库。但是，我不熟悉很难转换它的计算。所以，希望有一种方法可以通过一些修改来使用它。谢谢！

任何 cmets 都会有所帮助。

【问题讨论】：

简而言之，没有。原因：RcppNumerical 是面向 Eigen 库的。因此，MFuncGrad class 函数的基本定义是 Eigen::Ref<Eigen::VectorXd> 和 Eigen::Ref<const Eigen::VectorXd>

@coatless 感谢您的评论！现在，我要开始研究 Eigen 库了。

【参考方案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;
                    ^~~~~~

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

如果您使用相同的签名，事情应该会解决。所需的类型是根据 Eigen 对象定义的：

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

因此，您必须在 Aramdillo 和 Eigen 构造之间来回转换。

【讨论】：

感谢您的评论！正如你所说，我需要在 Aramdillo 和 Eigen 构造之间来回转换。

@WoojungBae 查看this answer 以获得有关转换的一些指示。这样做可能无需来回复制数据。不过，将您的代码转换为 Eigen 可能更有效。

谢谢！我希望不用转换上面所有的代码也可以。