2016-08-08 29 views
1

Tief im Inneren eines MCMC Algorithmus muss ich mit einem Vektor, also das folgende Stück RCPP und RcppArmadillo Code, um eine vom Benutzer bereitgestellte Liste von Matrizen multiplizieren wird mehrmals pro MCMC Iteration genannt:eine Liste von RcppArmadillo Matrizen

List mat_vec1 (const List& Mats, const vec& y) { 
    int n_list = Mats.size(); 
    Rcpp::List out(n_list); 
    for (int i = 0; i < n_list; ++i) { 
     out[i] = as<mat>(Mats[i]) * y; 
    } 
    return(out); 
} 

Die vom Benutzer bereitgestellte Liste Mats bleibt während der MCMC-, Vektor y Änderungen in jeder Iteration unverändert. Effizienz ist vorrangig und ich versuche zu sehen, ob ich den Code beschleunigen kann, indem ich die Elemente von Mats nicht oft zu arma :: mat konvertieren muss (es muss nur einmal gemacht werden). Ich habe versucht, den folgenden Ansatz

List arma_Mats (const List& Mats) { 
    int n_list = Mats.size(); 
    Rcpp::List res(n_list); 
    for (int i = 0; i < n_list; ++i) { 
     res[i] = as<mat>(Mats[i]); 
    } 
    return(res); 
} 

und dann

List mat_vec2 (const List& Mats, const vec& y) { 
    int n_list = Mats.size(); 
    Rcpp::List aMats = arma_Mats(Mats); 
    Rcpp::List out(n_list); 
    for (int i = 0; i < n_list; ++i) { 
     out[i] = aMats[i] * y; 
    } 
    return(out); 
} 

aber das scheint nicht zu funktionieren. Hinweise auf alternative/bessere Lösungen sind sehr willkommen.

+0

'List' ist ein R (kompatibler Rcpp) -Typ, der andere R (compatibale Rcpp) -Typen speichern kann (und solche, für die wir Konverter haben) - also müssen Sie' 'wrae()' Ihre 'arma: : Mat 'zuerst, um den Compiler mit einer 'SEXP'-kompatiblen Variable zu helfen. –

+0

Danke für den Zeiger. Wenn ich also richtig verstehe, wird 'Rcpp :: List' die Armadillo-Matrizen in einfache Rcpp-Matrizen umwandeln. Gibt es eine Möglichkeit, eine "Liste" von Armadillo-Matrizen zu erstellen, dann kann ich sie mit einem Armadillo-Vektor multiplizieren **, ohne dass ich sie wieder konvertieren muss? – user2692802

+0

Werfen Sie einen Blick auf die Armadillo-Dokumentation; 'Felder' Klassen können sein, was Sie wollen. Wenn Sie Daten an R zurücksenden möchten, benötigen Sie Rcpp-Typen oder Klassen, von denen R und Rcpp wissen. Dazu gehören Armadillo-Typen über RcppArmadillo. Studiere einige Beispiele ... –

Antwort

3

Ok, ich schrieb im Grunde die Antwort in dem Kommentar aber es kam dann zu mir, dass wir bereits ein funktionierendes Beispiel in den von RcppArmadillo.package.skeleton() erstellt Stub bieten:

// [[Rcpp::export]] 
Rcpp::List rcpparma_bothproducts(const arma::colvec & x) { 
    arma::mat op = x * x.t(); 
    double ip = arma::as_scalar(x.t() * x); 
    return Rcpp::List::create(Rcpp::Named("outer")=op, 
           Rcpp::Named("inner")=ip); 
} 

Dies ist eine Liste gibt das äußeree Produkt (a Matrix) und das innere Produkt (ein Skalar) des gegebenen Vektors.

Was schnell ist und was nicht: Ich empfehle, nicht zu vermuten, sondern Profil und messen Sie so viel wie Sie können. Meine Neigung wäre, mehr (eigenständigen) C++ - Code in Armadillo zu machen und nur am Ende zurückzukehren, um die Conversions zu minimieren.