C'è qualche funzione di ottimizzazione in Rcpp

Question

Quello che segue è il mio codice Rcpp e voglio ridurre al minimo la funzione obiettivo logtpoi (x, theta) rispetto a theta in R da 'nlminb'. Ho scoperto che è lento. Ho due domande:

1. Chiunque può migliorare il mio codice Rcpp? Grazie mille.
2. Esistono funzioni di ottimizzazione in Rcpp? Se sì, forse posso usarli direttamente in Rcpp. E come usarli? Grazie mille.

Il mio codice:

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

// [[Rcpp::depends("RcppArmadillo")]] 
// [[Rcpp::export]]
List dtpoi0(const IntegerVector& x, const NumericVector& theta){
//x is 3-dim vector; theta is a 6-dim parameter vector.
//be careful the order of theta1,...,theta6.
  double theta1 = theta[0]; double theta2 = theta[1];
  double theta3 = theta[2]; double theta4 = theta[3];
  double theta5 = theta[4]; double theta6 = theta[5];
  int x1 = x[0]; int x2 = x[1]; int x3 = x[2];
  IntegerVector z1 = IntegerVector::create(x1,x2); 
  IntegerVector z2 = IntegerVector::create(x1,x3);
  IntegerVector z3 = IntegerVector::create(x2,x3); 
  int s1 = min(z1); int s2 = min(z2); int s3 = min(z3);
  arma::imat missy(1,3,fill::zeros); arma::irowvec ijk={0,0,0};
  for (int i = 0; i <= s1; ++i) {
    for (int j = 0; j <= s2; ++j) {
      for (int k = 0; k <= s3; ++k) {
        if ((i+j <= s1) & (i+k <= s2) & ( j+k <= s3))
        { ijk = {i,j,k};
          missy = join_cols(missy,ijk);}
       }
     }
   }
  IntegerMatrix misy = as<IntegerMatrix>(wrap(missy));
  IntegerVector u1 = IntegerVector::create(0);
  IntegerVector u2 = IntegerVector::create(0);
  IntegerVector u3 = IntegerVector::create(0);
  IntegerVector u4 = IntegerVector::create(0);
  IntegerVector u5 = IntegerVector::create(0);
  IntegerVector u6 = IntegerVector::create(0);
  int total = misy.nrow();
  double fvalue = 0; 
  NumericVector part1(1); NumericVector part2(1); 
  NumericVector part3(1); NumericVector part4(1); 
  NumericVector part5(1); NumericVector part6(1);
  for (int l = 1; l < total; ++l) {
    u1 = IntegerVector::create(x1-misy(l,0)-misy(l,1));
    u2 = IntegerVector::create(x2-misy(l,0)-misy(l,2));
    u3 = IntegerVector::create(x3-misy(l,1)-misy(l,2));
    u4 = IntegerVector::create(misy(l,0));
    u5 = IntegerVector::create(misy(l,1));
    u6 = IntegerVector::create(misy(l,2));
    part1 = dpois(u1,theta1);
    part2 = dpois(u2,theta2);
    part3 = dpois(u3,theta3);
    part4 = dpois(u4,theta4);
    part5 = dpois(u5,theta5);
    part6 = dpois(u6,theta6);
    fvalue = fvalue + (part1*part2*part3*part4*part5*part6)[0]; }
   return(List::create(Named("misy") = misy,Named("fvalue") = fvalue)); 
  }
 
 // [[Rcpp::export]]
 NumericVector dtpoi(const IntegerMatrix& x, const NumericVector& theta){ 
//x is n*3 matrix, n is the number of observations.
   int n = x.nrow();
   NumericVector density(n);
   for (int i = 0; i < n; ++i){
     density(i) = dtpoi0(x.row(i),theta)["fvalue"];
   }
   return(density);
 }

// [[Rcpp::export]]
double logtpoi0(const IntegerMatrix& x,const NumericVector theta){
// theta must be a 6-dimiension parameter.
  double nln = -sum(log( dtpoi(x,theta) + 1e-60 ));
  if(arma::is_finite(nln)) {nln = nln;} else {nln = -1e10;} 
  return(nln); 
}

Answer 1

Enorme avvertimento in vista: non conosco davvero Armadillo. Ma ci ho provato perché il codice sembra interessante.

Alcune cose generali:

1. Non è necessario dichiarare le cose prima di assegnarle per la prima volta. In particolare, generalmente non è necessario dichiarare i vettori all'esterno di un ciclo se vengono utilizzati solo all'interno del ciclo. Questo probabilmente non è meno efficiente che dichiararli all'interno del ciclo. Tuttavia, se il codice è troppo lento, ha senso profilarlo attentamente e verificare se il presupposto è valido.
2. Molte delle tue dichiarazioni sono solo alias per elementi vettoriali e non sembrano necessarie.
3. I tuoi z{1…3}vettori non sono necessari. C ++ ha una minfunzione per trovare il minimo di due elementi.
4. dtpoi0contiene due loop principali. Entrambi sono stati pesantemente modificati nel mio codice:
   1. Il primo ciclo itera su molti messaggi kche non possono essere mai utilizzati, a causa dell'interno ifche verifica se i + jsupera s2. Inserendo questo controllo nella condizione di ciclo di j, eseguiamo meno kcicli.
      1. I tuoi ifusi &invece di &&. Come in R, usare &&piuttosto che &causare cortocircuiti . Anche se questo probabilmente non è più efficiente in questo caso, l'uso &&è idiomatico, mentre &provoca grattacapi (il mio codice usa andche è un modo alternativo di ortografia &&in C ++; preferisco la sua leggibilità).
   2. Il secondo loop esegue efficacemente un'operazione di matrice manualmente. Sento che dovrebbe esserci un modo per esprimerlo puramente con operazioni di matrice - ma come detto non sono un utente di Armadillo. Tuttavia, le mie modifiche tentano di vettorializzare il più possibile questa operazione (se non altro questo rende il codice più breve). Il dpoisprodotto interno è purtroppo ancora all'interno di un ciclo.
5. La logica di logtpoi0può essere resa più idiomatica e (IMHO) più leggibile utilizzando l'operatore condizionale invece di if.
6. la correttezza const è un grosso problema in C ++, poiché elimina le modifiche accidentali. Usalo constliberamente quando dichiari variabili che non dovrebbero cambiare.
7. In termini di efficienza, il problema maggiore quando si chiama dtpoio logtpoi0è probabilmente la conversione di missyin misy, che causa allocazioni e copie di memoria. Converti in solo IntegerMatrixquando necessario, ovvero quando restituisci effettivamente quel valore a R. Per questo motivo, ho diviso dtpoi0in due parti.
8. Un'altra inefficienza è data dal fatto che il primo ciclo in fa dtpoi0crescere una matrice aggiungendo colonne. Questo è un grande no-no. Tuttavia, riscrivere il codice per evitare ciò non è banale.

#include <algorithm>

#include <RcppArmadillo.h>

// [[Rcpp::depends("RcppArmadillo")]]

using namespace Rcpp;
using namespace arma;

imat dtpoi0_mat(const IntegerVector& x) {
    const int s1 = std::min(x[0], x[1]);
    const int s2 = std::min(x[0], x[2]);
    const int s3 = std::min(x[1], x[2]);
    imat missy(1, 3, fill::zeros);

    for (int i = 0; i <= s1; ++i) {
        for (int j = 0; j <= s2 and i + j <= s1; ++j) {
            for (int k = 0; k <= s3 and i + k <= s2 and j + k <= s3; ++k) {
                missy = join_cols(missy, irowvec{i, j, k});
            }
        }
    }

    return missy;
}

double dtpoi0_fvalue(const IntegerVector& x, const NumericVector& theta, imat& missy) {
    double fvalue = 0.0;
    ivec xx = as<ivec>(x);
    missy.each_row([&](irowvec& v) {
        const ivec u(join_cols(xx - v(uvec{0, 0, 1}) - v(uvec{1, 2, 3}), v));
        double prod = 1;
        for (int i = 0; i < u.n_elem; ++i) {
            prod *= R::dpois(u[i], theta[i], 0);
        }
        fvalue += prod;
    });
    return fvalue;
}

double dtpoi0_fvalue(const IntegerVector& x, const NumericVector& theta) {
    imat missy = dtpoi0_mat(x);
    return dtpoi0_fvalue(x, theta, missy);
}

// [[Rcpp::export]]
List dtpoi0(const IntegerVector& x, const NumericVector& theta) {
    imat missy = dtpoi0_mat(x);
    const double fvalue = dtpoi0_fvalue(x, theta, missy);
    return List::create(Named("misy") = as<IntegerMatrix>(wrap(missy)), Named("fvalue") = fvalue);
}

// [[Rcpp::export]]
NumericVector dtpoi(const IntegerMatrix& x, const NumericVector& theta) {
    //x is n*3 matrix, n is the number of observations.
    int n = x.nrow();
    NumericVector density(n);

    for (int i = 0; i < n; ++i){
        density(i) = dtpoi0_fvalue(x.row(i), theta);
    }

    return density;
}

// [[Rcpp::export]]
double logtpoi0(const IntegerMatrix& x, const NumericVector theta) {
    // theta must be a 6-dimension parameter.
    const double nln = -sum(log(dtpoi(x, theta) + 1e-60));
    return is_finite(nln) ? nln : -1e10;
}

Importante: viene compilato, ma non posso testarne la correttezza. È del tutto possibile (anzi probabile!) Che il mio refactor abbia introdotto errori. Dovrebbe quindi essere visto solo come uno schizzo della soluzione e non dovrebbe in alcun modo essere copiato e incollato in un'applicazione.