optimization - 带有 Rcpparmadillo 的 Rcppnumerical 和 RcppEigen 中的 obtim_lbfgs 函数
问题描述
[编辑]
我正在寻找在 Rcppnumerical 和 RcppEigen 中使用 optim_lbfgs 函数和 Rcpparmadillo 的方法。我遵循Rcpp Integration for Numerical Computing Libraries中的方法,但它没有解决错误,cannot declare variable 'obj' to be of abstract type 'scoreftn_mns'
. 但是现在,我修复了一些代码以使其正常工作。我在Rcppnumericalbeta
中定义as并在模板(?)中将其转换为 arma::vec。Eigen::VectorXd beta(p);
这是我正在尝试做的代码。
// [[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 库。但是,我不熟悉很难转换它的计算。所以,希望有一种方法可以通过一些修改来使用它。谢谢!
任何评论都会有所帮助。
解决方案
我在编译您的代码时收到的其余错误消息很有趣:
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 构造之间来回转换。
推荐阅读
- reactjs - Next.js - 关于 dedupeInterval 和 refreshInterval 的 SWR 钩子问题
- java - 以堆栈的相反顺序显示字母
- python - 图例句柄在 matplotlib 中未正确对齐
- python - 使用过去三个月的第一天和最后一天生成数据框
- r - 将旧的、奇怪的格式、文本文件导入 R
- bloomberg - 彭博商品期限结构数据
- javascript - 如何在 React 中动态导入组件?
- python - MS Azure Web 应用程序上的 Django-Python 日志错误
- google-cloud-platform - 无法删除 GCP 文件夹,因为它包含活动资源
- python - Python 响应侦听器,直到响应包含关键字/值