Saat ini saya sedang mengerjakan metode Bayesian yang membutuhkan beberapa langkah optimasi model multinomial logit per iterasi. Saya menggunakan optim () untuk melakukan optimasi tersebut, dan fungsi objektif yang ditulis dalam profil R. mengungkapkan bahwa optim () adalah hambatan utama.

Setelah menggali sekitar, saya menemukan pertanyaan ini di mana mereka menyarankan bahwa pengodean ulang fungsi tujuan dengan Rcppdapat mempercepat proses. Saya mengikuti saran itu dan mencatat kembali fungsi obyektif saya Rcpp, tetapi akhirnya menjadi lebih lambat (sekitar dua kali lebih lambat!).

Ini adalah pertama kalinya saya dengan Rcpp(atau apa pun yang berhubungan dengan C ++) dan saya tidak dapat menemukan cara vektorisasi kode. Adakah yang tahu cara membuatnya lebih cepat?

Tl; dr: Implementasi fungsi saat ini di Rcpp tidak secepat R vektor; bagaimana membuatnya lebih cepat?

Contoh yang dapat direproduksi :

1) Menentukan fungsi obyektif dalam Rdan Rcpp: log-kemungkinan model intersep hanya multinomial

library(Rcpp)
library(microbenchmark)

llmnl_int <- function(beta, Obs, n_cat) {
  n_Obs     <- length(Obs)
  Xint      <- matrix(c(0, beta), byrow = T, ncol = n_cat, nrow = n_Obs)
  ind       <- cbind(c(1:n_Obs), Obs)
  Xby       <- Xint[ind]
  Xint      <- exp(Xint)
  iota      <- c(rep(1, (n_cat)))
  denom     <- log(Xint %*% iota)
  return(sum(Xby - denom))
}

cppFunction('double llmnl_int_C(NumericVector beta, NumericVector Obs, int n_cat) {

    int n_Obs = Obs.size();

    NumericVector betas = (beta.size()+1);
    for (int i = 1; i < n_cat; i++) {
        betas[i] = beta[i-1];
    };

    NumericVector Xby = (n_Obs);
    NumericMatrix Xint(n_Obs, n_cat);
    NumericVector denom = (n_Obs);
    for (int i = 0; i < Xby.size(); i++) {
        Xint(i,_) = betas;
        Xby[i] = Xint(i,Obs[i]-1.0);
        Xint(i,_) = exp(Xint(i,_));
        denom[i] = log(sum(Xint(i,_)));
    };

    return sum(Xby - denom);
}')

2) Bandingkan efisiensinya:

## Draw sample from a multinomial distribution
set.seed(2020)
mnl_sample <- t(rmultinom(n = 1000,size = 1,prob = c(0.3, 0.4, 0.2, 0.1)))
mnl_sample <- apply(mnl_sample,1,function(r) which(r == 1))

## Benchmarking
microbenchmark("llmml_int" = llmnl_int(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
               "llmml_int_C" = llmnl_int_C(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
               times = 100)
## Results
# Unit: microseconds
#         expr     min       lq     mean   median       uq     max neval
#    llmnl_int  76.809  78.6615  81.9677  79.7485  82.8495 124.295   100
#  llmnl_int_C 155.405 157.7790 161.7677 159.2200 161.5805 201.655   100

3) Sekarang memanggil mereka optim:

## Benchmarking with optim
microbenchmark("llmnl_int" = optim(c(4,2,1), llmnl_int, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
               "llmnl_int_C" = optim(c(4,2,1), llmnl_int_C, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
               times = 100)
## Results
# Unit: milliseconds
#         expr      min       lq     mean   median       uq      max neval
#    llmnl_int 12.49163 13.26338 15.74517 14.12413 18.35461 26.58235   100
#  llmnl_int_C 25.57419 25.97413 28.05984 26.34231 30.44012 37.13442   100

Saya agak terkejut bahwa implementasi vectorised di R lebih cepat. Menerapkan versi yang lebih efisien dalam Rcpp (katakanlah, dengan RcppArmadillo?) Dapat menghasilkan keuntungan? Apakah ide yang lebih baik untuk mengode ulang semua yang ada di Rcpp menggunakan pengoptimal C ++?

PS: posting pertama kali di Stackoverflow!

Secara umum jika Anda dapat menggunakan fungsi vektor, Anda akan mendapati (hampir) secepat menjalankan kode Anda langsung di Rcpp. Ini karena banyak fungsi vektor dalam R (hampir semua fungsi vektor dalam Basis R) ditulis dalam C, Cpp atau Fortran dan karena itu sering ada sedikit untuk mendapatkan.

Yang mengatakan, ada perbaikan untuk mendapatkan kode Anda Rdan Anda Rcpp. Optimalisasi berasal dari mempelajari kode dengan cermat, dan menghapus langkah yang tidak perlu (penugasan memori, jumlah, dll.).

Mari kita mulai dengan Rcppoptimasi kode.

Dalam kasus Anda, optimasi utama adalah untuk menghapus perhitungan matriks dan vektor yang tidak perlu. Kode ini pada intinya

1. Pergeseran beta
2. menghitung log dari jumlah exp (shift beta) [log-sum-exp]
3. gunakan Obs sebagai indeks untuk beta yang dipindahkan dan jumlahkan semua probabilitas
4. kurangi log-sum-exp

Dengan menggunakan pengamatan ini, kami dapat mengurangi kode Anda menjadi 2 for-loop. Perhatikan bahwa sumini hanyalah for-loop (lebih atau kurang for(i = 0; i < max; i++){ sum += x }:) sehingga menghindari penjumlahan dapat mempercepat kode yang lebih jauh (dalam kebanyakan situasi ini adalah optimasi yang tidak perlu!). Selain itu, input Anda Obsadalah vektor integer, dan kami dapat lebih mengoptimalkan kode dengan menggunakan IntegerVectortipe untuk menghindari penuangan doubleelemen ke integernilai (kredit ke jawaban Ralf Stubner).

cppFunction('double llmnl_int_C_v2(NumericVector beta, IntegerVector Obs, int n_cat)
 {

    int n_Obs = Obs.size();

    NumericVector betas = (beta.size()+1);
    //1: shift beta
    for (int i = 1; i < n_cat; i++) {
        betas[i] = beta[i-1];
    };
    //2: Calculate log sum only once:
    double expBetas_log_sum = log(sum(exp(betas)));
    // pre allocate sum
    double ll_sum = 0;

    //3: Use n_Obs, to avoid calling Xby.size() every time 
    for (int i = 0; i < n_Obs; i++) {
        ll_sum += betas(Obs[i] - 1.0) ;
    };
    //4: Use that we know denom is the same for all I:
    ll_sum = ll_sum - expBetas_log_sum * n_Obs;
    return ll_sum;
}')

Perhatikan bahwa saya telah menghapus beberapa alokasi memori dan menghapus perhitungan yang tidak perlu di for-loop. Juga saya telah menggunakan yang denomsama untuk semua iterasi dan hanya dikalikan untuk hasil akhir.

Kami dapat melakukan optimasi serupa di R-code Anda, yang menghasilkan fungsi di bawah ini:

llmnl_int_R_v2 <- function(beta, Obs, n_cat) {
    n_Obs <- length(Obs)
    betas <- c(0, beta)
    #note: denom = log(sum(exp(betas)))
    sum(betas[Obs]) - log(sum(exp(betas))) * n_Obs
}

Perhatikan kompleksitas fungsi telah berkurang secara drastis sehingga lebih mudah dibaca orang lain. Hanya untuk memastikan bahwa saya belum mengacaukan kode di suatu tempat mari kita periksa apakah mereka mengembalikan hasil yang sama:

set.seed(2020)
mnl_sample <- t(rmultinom(n = 1000,size = 1,prob = c(0.3, 0.4, 0.2, 0.1)))
mnl_sample <- apply(mnl_sample,1,function(r) which(r == 1))

beta = c(4,2,1)
Obs = mnl_sample 
n_cat = 4
xr <- llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xr2 <- llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc <- llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc2 <- llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
all.equal(c(xr, xr2), c(xc, xc2))
TRUE

baik itu melegakan.

Kinerja:

Saya akan menggunakan microbenchmark untuk menggambarkan kinerja. Fungsi yang dioptimalkan cepat, jadi saya akan menjalankan fungsi 1e5kali untuk mengurangi efek pengumpul sampah

microbenchmark("llmml_int_R" = llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmml_int_C" = llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmnl_int_R_v2" = llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmml_int_C_v2" = llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               times = 1e5)
#Output:
#Unit: microseconds
#           expr     min      lq       mean  median      uq        max neval
#    llmml_int_R 202.701 206.801 288.219673 227.601 334.301  57368.902 1e+05
#    llmml_int_C 250.101 252.802 342.190342 272.001 399.251 112459.601 1e+05
# llmnl_int_R_v2   4.800   5.601   8.930027   6.401   9.702   5232.001 1e+05
# llmml_int_C_v2   5.100   5.801   8.834646   6.700  10.101   7154.901 1e+05

Di sini kita melihat hasil yang sama seperti sebelumnya. Sekarang fungsi-fungsi baru kira-kira 35x lebih cepat (R) dan 40x lebih cepat (Cpp) dibandingkan dengan bagian-bagian pertama. Yang cukup menarik, Rfungsi yang dioptimalkan masih sangat sedikit (0,3 ms atau 4%) lebih cepat daripada Cppfungsi yang dioptimalkan . Taruhan terbaik saya di sini adalah bahwa ada beberapa overhead dari Rcpppaket, dan jika ini dihapus keduanya akan identik atau R.

Demikian pula kita dapat memeriksa kinerja menggunakan Optim.

microbenchmark("llmnl_int" = optim(beta, llmnl_int, Obs = mnl_sample, 
                                   n_cat = n_cat, method = "BFGS", hessian = F, 
                                   control = list(fnscale = -1)),
               "llmnl_int_C" = optim(beta, llmnl_int_C, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               "llmnl_int_R_v2" = optim(beta, llmnl_int_R_v2, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               "llmnl_int_C_v2" = optim(beta, llmnl_int_C_v2, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               times = 1e3)
#Output:
#Unit: microseconds
#           expr       min        lq      mean    median         uq      max neval
#      llmnl_int 29541.301 53156.801 70304.446 76753.851  83528.101 196415.5  1000
#    llmnl_int_C 36879.501 59981.901 83134.218 92419.551 100208.451 190099.1  1000
# llmnl_int_R_v2   667.802  1253.452  1962.875  1585.101   1984.151  22718.3  1000
# llmnl_int_C_v2   704.401  1248.200  1983.247  1671.151   2033.401  11540.3  1000

Sekali lagi hasilnya sama.

Kesimpulan:

Sebagai kesimpulan singkat perlu dicatat bahwa ini adalah salah satu contoh, di mana mengkonversi kode Anda ke Rcpp tidak terlalu sepadan dengan masalahnya. Ini tidak selalu terjadi, tetapi sering kali perlu melihat fungsi Anda, untuk melihat apakah ada area kode Anda, di mana perhitungan yang tidak perlu dilakukan. Terutama dalam situasi di mana seseorang menggunakan fungsi vektor buildin, seringkali tidak sepadan dengan waktu untuk mengkonversi kode ke Rcpp. Lebih sering orang dapat melihat peningkatan besar jika menggunakan for-loopskode yang tidak dapat dengan mudah di-vektor-kan untuk menghapus for-loop.

Fungsi C ++ Anda dapat dibuat lebih cepat menggunakan pengamatan berikut. Setidaknya yang pertama juga dapat digunakan dengan fungsi R Anda:

• Cara Anda menghitung denom[i]sama untuk setiap i. Karena itu masuk akal untuk menggunakan double denomdan melakukan perhitungan ini hanya sekali. Saya juga faktor keluar mengurangi istilah umum ini pada akhirnya.

• Pengamatan Anda sebenarnya adalah vektor integer di sisi R, dan Anda menggunakannya sebagai integer di C ++ juga. Menggunakan a IntegerVectoruntuk memulai dengan membuat banyak casting tidak perlu.

• Anda dapat mengindeks NumericVectormenggunakan IntegerVectorC ++ juga. Saya tidak yakin apakah ini membantu kinerja, tetapi membuat kode sedikit lebih pendek.

• Beberapa perubahan yang lebih terkait dengan gaya daripada kinerja.

Hasil:

double llmnl_int_C(NumericVector beta, IntegerVector Obs, int n_cat) {

    int n_Obs = Obs.size();

    NumericVector betas(beta.size()+1);
    for (int i = 1; i < n_cat; ++i) {
        betas[i] = beta[i-1];
    };

    double denom = log(sum(exp(betas)));
    NumericVector Xby = betas[Obs - 1];

    return sum(Xby) - n_Obs * denom;
}

Bagi saya fungsi ini kira-kira sepuluh kali lebih cepat daripada fungsi R Anda.

Saya dapat memikirkan empat optimisasi potensial atas jawaban Ralf dan Olivers.

(Anda harus menerima jawaban mereka, tetapi saya hanya ingin menambahkan 2 sen saya).

1) Gunakan // [[Rcpp::export(rng = false)]]sebagai tajuk komentar ke fungsi dalam file C ++ terpisah. Ini menyebabkan kecepatan ~ 80% pada mesin saya. (Ini adalah saran paling penting dari 4).

2) Lebih suka cmathbila memungkinkan. (Dalam hal ini, sepertinya tidak ada bedanya).

3) Hindari alokasi kapan pun memungkinkan, misalnya jangan beralih betake vektor baru.

4) Regangkan tujuan: gunakan SEXPparameter daripada vektor Rcpp. (Ditinggalkan sebagai latihan untuk pembaca). Vektor Rcpp adalah pembungkus yang sangat tipis, tetapi mereka masih pembungkus dan ada overhead kecil.

Saran-saran ini tidak akan penting, jika bukan karena Anda memanggil fungsi dalam loop ketat optim. Jadi setiap overhead sangat penting.

Bangku:

microbenchmark("llmnl_int_R_v1" = optim(beta, llmnl_int, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_R_v2" = optim(beta, llmnl_int_R_v2, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v2" = optim(beta, llmnl_int_C_v2, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v3" = optim(beta, llmnl_int_C_v3, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v4" = optim(beta, llmnl_int_C_v4, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             times = 1000)


Unit: microseconds
expr      min         lq       mean     median         uq        max neval cld
llmnl_int_R_v1 9480.780 10662.3530 14126.6399 11359.8460 18505.6280 146823.430  1000   c
llmnl_int_R_v2  697.276   735.7735  1015.8217   768.5735   810.6235  11095.924  1000  b 
llmnl_int_C_v2  997.828  1021.4720  1106.0968  1031.7905  1078.2835  11222.803  1000  b 
llmnl_int_C_v3  284.519   295.7825   328.5890   304.0325   328.2015   9647.417  1000 a  
llmnl_int_C_v4  245.650   256.9760   283.9071   266.3985   299.2090   1156.448  1000 a 

v3 adalah jawaban Oliver rng=false. v4 dengan Saran # 2 dan # 3 disertakan.

Fungsi:

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;

// [[Rcpp::export(rng = false)]]
double llmnl_int_C_v4(NumericVector beta, IntegerVector Obs, int n_cat) {

  int n_Obs = Obs.size();
  //2: Calculate log sum only once:
  // double expBetas_log_sum = log(sum(exp(betas)));
  double expBetas_log_sum = 1.0; // std::exp(0)
  for (int i = 1; i < n_cat; i++) {
    expBetas_log_sum += std::exp(beta[i-1]);
  };
  expBetas_log_sum = std::log(expBetas_log_sum);

  double ll_sum = 0;
  //3: Use n_Obs, to avoid calling Xby.size() every time 
  for (int i = 0; i < n_Obs; i++) {
    if(Obs[i] == 1L) continue;
    ll_sum += beta[Obs[i]-2L];
  };
  //4: Use that we know denom is the same for all I:
  ll_sum = ll_sum - expBetas_log_sum * n_Obs;
  return ll_sum;
}