Эффективность и скорость R-кода с использованием Rcpp

Question

Это сообщение об ускорении R-кода с использованием пакета Rcpp, чтобы избежать рекурсивных циклов.

Мой вход определяют на следующем примере (длина 7), которая является частью data.frame (длина 51673), который я использовал:

S=c(906.65,906.65,906.65,906.65,906.65,906.65,906.65) 
T=c(0.1371253,0.1457896,0.1248953,0.1261278,0.1156931,0.0985253,0.1332596) 
r=c(0.013975,0.013975,0.013975,0.013975,0.013975,0.013975,0.013975)    
h=c(0.001332596,0.001248470,0.001251458,0.001242143,0.001257921,0.0,0.0)  
P=c(3,1,5,2,1,4,2) 
A= data.frame(S=S,T=T,r=r,h=h,P=P) 

     S   T  r   h Per 
1 906.65 0.1971253 0.013975 0.001332596 3 
2 906.65 0.1971253 0.013975 0.001248470 1 
3 906.65 0.1971253 0.013975 0.001251458 5 
4 906.65 0.1971253 0.013975 0.001242143 2 
5 906.65 0.1971253 0.013975 0.001257921 1 
6 906.65 0.1971253 0.013975 0.0
7 906.65 0.1971253 0.013975 0.0

Параметры:

w=0.001; b=0.2; a=0.0154; c=0.0000052; neta=-0.70 

У меня есть следующий код функции, которую я хочу использовать:

F<-function(x,w,b,a,c,neta,S,T,r,P){ 
    u=1i*x 
    nu=(1/(neta^2))*(((1-2*neta)^(1/2))-1) 
    # Recursion back to time t 
    # Terminal condition for the A and B 
    A_Q=0 
    B_Q=0 
    steps<-round(T*250,0) 

    for (j in 1:steps){ 
     A_Q= A_Q+ r*u + w*B_Q-(1/2)*log(1-2*a*(neta^4)*B_Q) 
     B_Q= b*B_Q+u*nu+ (1/neta^2)*(1-sqrt((1-2*a*(neta^4)*B_Q)*(1- 2*c*B_Q - 2*u*neta))) 
    } 
    F= exp(log(S)*u + A_Q + B_Q*h[P]) 
    return(F) 
} 

S = A$S ; r= A$r ; T= A$T ; P=A$P; h= A$h 

Затем я хочу применить предыдущий АЯ функция с помощью моего Data.set вектора длины N = 100000:

Z=length(S); N=100000 ; alpha=2 ; delta= 0.25  
    lambda=(2*pi)/(N*delta) 

res = matrix(nrow=N, ncol=Z) 
    for (i in 1:N){ 
    for (j in 1:Z){ 
     res[i,j]= Re(F(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j])) 
    } 
    } 

Но это занимает много времени: она занимает 20 секунд, чтобы выполнить эту строку коды для N = 100, но я хочу выполните его для N = 100000 раз, общее время выполнения может занять несколько часов. Как точно настроить вышеуказанный код с помощью Rcpp, чтобы сократить время выполнения и получить эффективную программу?

Возможно ли сократить время выполнения, и если да, предложите мне решение даже без Rcpp.

Спасибо.

Answer 1

Ваша функция F может быть преобразован в C++ довольно легко воспользовавшись vec и cx_vec классов в Armadillo library (доступ через RcppArmadillo пакет) - который имеет большую поддержку векторизованных расчетов.

---

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)]] 

// [[Rcpp::export]] 
arma::cx_vec Fcpp(const arma::cx_vec& x, double w, double b, double a, double c, 
        double neta, const arma::vec& S, const arma::vec& T, 
        const arma::vec& r, Rcpp::IntegerVector P, Rcpp::NumericVector h) { 

    arma::cx_vec u = x * arma::cx_double(0.0,1.0); 

    double nu = (1.0/std::pow(neta,2.0)) * (std::sqrt(1.0-2.0*neta)-1.0); 
    arma::cx_vec A_Q(r.size()); 
    arma::cx_vec B_Q(r.size()); 
    arma::vec steps = arma::round(T*250.0); 

    for (size_t j = 0; j < steps.size(); j++) { 
    for (size_t k = 0; k < steps[j]; k++) { 
     A_Q = A_Q + r*u + w*B_Q - 
       0.5*arma::log(1.0 - 2.0*a*std::pow(neta,4.0)*B_Q); 
     B_Q = b*B_Q + u*nu + (1.0/std::pow(neta,2.0)) * 
       (1.0 - arma::sqrt((1.0 - 2.0*a*std::pow(neta,4.0)*B_Q) * 
       (1.0 - 2.0*c*B_Q - 2.0*u*neta))); 
    } 
    } 

    arma::vec hP = Rcpp::as<arma::vec>(h[P-1]); 
    arma::cx_vec F = arma::exp(arma::log(S)*u + A_Q + B_Q*hP); 

    return F; 
} 

---

Только несколько незначительных изменений отметить:

• Я использую arma:: функции для векторизованных расчетов, таких как arma::log, arma::exp, arma::round, arma::sqrt, а также различные перегруженные операторы (*, +, -); но используя std::pow и std::sqrt для скалярных вычислений. В R это абстрагируется от нас, но здесь мы должны различать две ситуации.
• Ваша функция F имеет один цикл - for (i in 1:steps) - но версия C++ имеет два, только из-за различий в семантике цикла между двумя языками.
• Большинство из входных векторов arma:: классов (в отличие от использования Rcpp::NumericVector и Rcpp::ComplexVector), исключение составляет P и h, так как Rcpp векторы предлагают R-подобный доступ элемент - например, h[P-1]. Также обратите внимание, что P должно быть смещено на 1 (индексирование на основе 0 на C++), а затем преобразовано в вектор Armadillo (hP) с использованием Rcpp::as<arma::vec>, так как ваш компилятор будет жаловаться, если вы попытаетесь умножить cx_vec на NumericVector (B_Q*hP) ,
• Я добавил параметр функции h - не стоит полагаться на существование глобальной переменной h, которую вы делали в F. Если вам нужно использовать его в теле функции, вы должны передать его в функцию.

---

Я изменил название своей функции Fr и сделать бенчмаркинг немного проще, я просто завернутые свой двойной цикл, который заполнит матрицу res в функции Fr и Fcpp:

loop_Fr <- function(mat = res) { 
    for (i in 1:N) { 
    for (j in 1:Z) { 
     mat[i,j]= Re(Fr(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j],h)) 
    } 
    } 
    return(mat) 
} 
loop_Fcpp <- function(mat = res) { 
    for (i in 1:N) { 
    for (j in 1:Z) { 
     mat[i,j]= Re(Fcpp(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j],h)) 
    } 
    } 
    return(mat) 
} 
## 
R> all.equal(loop_Fr(),loop_Fcpp()) 
[1] TRUE 

---

Я сравнил две функции N = 100, N = 1000 и N = 100000 (который взял навсегда) - adjusti ng lambda и res соответственно, но сохраняя все остальное одинаково. Вообще говоря, Fcpp примерно в 10 раз быстрее, чем Fr на моем компьютере:

N <- 100 
lambda <- (2*pi)/(N*delta) 
res <- matrix(nrow=N, ncol=Z) 
## 
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=50L) 
Unit: milliseconds 
     expr  min  lq median  uq  max neval 
    loop_Fr() 142.44694 146.62848 148.97571 151.86318 186.67296 50 
loop_Fcpp() 14.72357 15.26384 15.58604 15.85076 20.19576 50 

---

N <- 1000 
lambda <- (2*pi)/(N*delta) 
res <- matrix(nrow=N, ncol=Z) 
## 
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=50L) 
Unit: milliseconds 
     expr  min  lq median  uq  max neval 
    loop_Fr() 1440.8277 1472.4429 1491.5577 1512.5636 1565.6914 50 
loop_Fcpp() 150.6538 153.2687 155.4156 158.0857 181.8452 50 

---

N <- 100000 
lambda <- (2*pi)/(N*delta) 
res <- matrix(nrow=N, ncol=Z) 
## 
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=2L) 
Unit: seconds 
     expr  min  lq median  uq  max neval 
    loop_Fr() 150.14978 150.14978 150.33752 150.52526 150.52526  2 
loop_Fcpp() 15.49946 15.49946 15.75321 16.00696 16.00696  2 

---

Другие переменные, представленные в вашем вопросе:

S <- c(906.65,906.65,906.65,906.65,906.65,906.65,906.65) 
T <- c(0.1371253,0.1457896,0.1248953,0.1261278,0.1156931,0.0985253,0.1332596) 
r <- c(0.013975,0.013975,0.013975,0.013975,0.013975,0.013975,0.013975)    
h <- c(0.001332596,0.001248470,0.001251458,0.001242143,0.001257921,0.0,0.0)  
P <- c(3,1,5,2,1,4,2) 
w <- 0.001; b <- 0.2; a <- 0.0154; c <- 0.0000052; neta <- (-0.70) 
Z <- length(S) 
alpha <- 2; delta <- 0.25