Estoy tratando de devolver un montón de matrices usando RCPP. Mi código a continuación es extremadamente ineficiente. Me gustaría saber si el siguiente código puede ser eficiente.

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

// [[Rcpp::export]]
Rcpp::List hello( 
    const arma::rowvec& g,
    const int& n, 
    const int& p,
    const arma::mat& S,
    const arma::mat& zc,
    const arma::rowvec& dl){
  Rcpp::List ht(n);

  for(int t=0; t < n;++t){

    arma::mat hhat(p,n);
    hhat.fill(0.0);
    for(int i = 0;i < n; ++i){
      arma::mat h(p,1);
      h.fill(0.0);
      if (t > i){
        for(int u=i;u <= t; ++u){
          arma::rowvec zr = zc.rows(i,i);
          h += exp(arma::as_scalar(g*zr.t())) * (zr.t() - S.cols(u,u))*dl(u);
        }
      }
      hhat.cols(i,i) = h;
    }
    ht[t] = hhat;
  }

  // Specify list length
  Rcpp::List res(1);
  res[0] = ht;

  return(res);
}

Aquí está el ejemplo.

g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(4800),nrow=3,ncol=1600)
dl=runif(1600)
z=matrix(runif(4800),nrow=1600,ncol=3)
ptm=proc.time();kkk= hello(g=g,n=n,p=p,S = S,zc=z,dl = dl);proc.time()-ptm;
 user  system elapsed 
  31.25    0.00   31.30 

Cualquier ayuda sería apreciada.

Siguiendo el código actualizado. Inicialmente estaba regresando la lista de una lista. Ahora devuelve una lista. Esto reduce el tiempo de computación en 10 segundos. Espero que este código se pueda mejorar aún más.

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

// [[Rcpp::export]]
Rcpp::List hello( 
    const arma::rowvec& g,
    const int& n, 
    const int& p,
    const arma::mat& S,
    const arma::mat& zc,
    const arma::rowvec& dl){
  Rcpp::List ht(n);


  for(int t=0; t < n;++t){

    arma::mat hhat(p,n);
    hhat.zeros();
    for(int i = 0;i < n; ++i){
      arma::mat h(p,1);
      // h.fill(0.0);
      h.zeros();
      if (t > i){
        for(int u=i;u <= t; ++u){
          //arma::rowvec zr = zc.rows(i,i);
          h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
        }
      }
      hhat.col(i) = h;
    }
    ht[t] = hhat;
  }

  // Specify list length
  // Rcpp::List res(1);
  // res[0] = ht;

  return(ht);
}

La fórmula que estoy tratando de implementar se proporciona a continuación. ingrese la descripción de la imagen aquí

1
Hello 20 dic. 2019 a las 01:49

2 respuestas

La mejor respuesta

En mi otra respuesta, miré la eficiencia de la devolución de datos y las optimizaciones simples. Aquí quiero ver algo diferente: la optimización del algoritmo.

Desea calcular hhat(i, t) para 0 <= i, t < n y i < t. Al observar su fórmula, vemos que la dependencia de hhat en i y t es muy diferente. En particular, hhat(i, t + 1) se puede escribir como hhat(i, t) + something. En este momento su ciclo externo ha terminado t y está volviendo a calcular todos estos valores intermedios. Al cambiar el orden de los bucles, es fácil realizar cada uno de estos cálculos una sola vez, reduciendo el algoritmo a dos bucles anidados. Esto significa que debe generar las matrices resultantes por separado. Y como no puede almacenar un arma::mat dentro de un Rcpp::List, necesito un std::vector adicional para el almacenamiento:

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

// [[Rcpp::export]]
Rcpp::List hello( 
        const arma::rowvec& g,
        const int& n, 
        const int& p,
        const arma::mat& S,
        const arma::mat& zc,
        const arma::rowvec& dl){

    std::vector<arma::mat> foo(n);
    for(int t=0; t < n;++t){
        arma::mat hhat(p,n);
        hhat.zeros();
        foo[t] = hhat;
    }

    for(int i = 0;i < n; ++i){
        arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
        for(int t=i+1; t < n;++t){
            h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
            foo[t].col(i) = h;
        }
    }

    Rcpp::List ht(n);
    for(int t=0; t < n;++t){
        ht[t] = foo[t];
    }

    return(ht);
}

// [[Rcpp::export]]
Rcpp::List hello_orig( 
        const arma::rowvec& g,
        const int& n, 
        const int& p,
        const arma::mat& S,
        const arma::mat& zc,
        const arma::rowvec& dl){
    Rcpp::List ht(n);


    for(int t=0; t < n;++t){

        arma::mat hhat(p,n);
        hhat.zeros();
        for(int i = 0;i < n; ++i){
            arma::mat h(p,1);
            h.zeros();
            if (t > i){
                for(int u=i;u <= t; ++u){
                    h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
                }
            }
            hhat.col(i) = h;
        }
        ht[t] = hhat;
    }

    return(ht);
}

/***R
g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(p*n),nrow=p,ncol=n)
dl=runif(n)
z=matrix(runif(p*n),nrow=n,ncol=p)
bench::mark(hello_orig(g=g,n=n,p=p,S = S,zc=z,dl = dl),
            hello(g=g,n=n,p=p,S = S,zc=z,dl = dl))
*/

Resultado:

# A tibble: 2 x 13
  expression                                                 min median `itr/sec` mem_alloc
  <bch:expr>                                              <bch:> <bch:>     <dbl> <bch:byt>
1 hello_orig(g = g, n = n, p = p, S = S, zc = z, dl = dl)  14.2s  14.2s    0.0703    58.7MB
2 hello(g = g, n = n, p = p, S = S, zc = z, dl = dl)      53.9ms 85.9ms   11.1       58.7MB
# … with 8 more variables: `gc/sec` <dbl>, n_itr <int>, n_gc <dbl>, total_time <bch:tm>,
#   result <list>, memory <list>, time <list>, gc <list>

¡Más de un factor 100 más rápido!

Puede obtener un código más limpio (y tal vez incluso un código un poco más rápido) haciendo sugerencias de @coatless en los comentarios para usar un arma::cube. Sin embargo, la forma más compacta le dará una estructura de retorno diferente. En lugar de una lista de p x n, obtendrá una matriz p x n x n:

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

// [[Rcpp::export]]
arma::cube coatless( 
        const arma::rowvec& g,
        const int& n, 
        const int& p,
        const arma::mat& S,
        const arma::mat& zc,
        const arma::rowvec& dl){

    arma::cube ht(p, n, n);
    ht.zeros();

    for(int i = 0;i < n; ++i){
        arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
        for(int t=i+1; t < n;++t){
            h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
            ht.slice(t).col(i) = h;
        }
    }

    return(ht);
}
3
Ralf Stubner 25 dic. 2019 a las 10:33

El título de su pregunta le hace pensar que ve el problema al devolver los datos a R. Tenga la seguridad de que esto no es un problema. Puede verificar esto fácilmente llamando a una función que devuelve matrices de ceros en el tamaño requerido:

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

// [[Rcpp::export]]
Rcpp::List minimal( 
        const arma::rowvec& g,
        const int& n, 
        const int& p,
        const arma::mat& S,
        const arma::mat& zc,
        const arma::rowvec& dl){
    Rcpp::List ht(n);

    for(int t=0; t < n;++t){

        arma::mat hhat(p,n);
        hhat.zeros();
        ht[t] = hhat;
    }

    return(ht);
}

En mi sistema, esta función tarda aproximadamente 0.01 s con sus datos de entrada. En otras palabras, su función real pasa la mayor parte del tiempo calculando los resultados reales.

En cuanto a la optimización de esa parte, sería útil si pudiera proporcionar una idea de lo que está tratando de implementar, p. con la ayuda de fórmulas matemáticas. Tal como está, solo puedo hacer algunos cambios simples:

  • En el ciclo i solo haces algo por t > i. Por lo tanto, es suficiente dejar que el ciclo se ejecute hasta i < t.

  • El bucle u puede formularse como un producto matriz-vector, para el cual existen implementaciones eficientes.

Con cambios como este termino con

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

// [[Rcpp::export]]
Rcpp::List hello( 
        const arma::rowvec& g,
        const int& n, 
        const int& p,
        const arma::mat& S,
        const arma::mat& zc,
        const arma::rowvec& dl){
    Rcpp::List ht(n);

    for(int t=0; t < n;++t){

        arma::mat hhat(p,n);
        hhat.zeros();
        for(int i = 0;i < t; ++i){
            arma::mat Sit = S.cols(i,t);
            hhat.col(i) = - exp(arma::as_scalar(g*zc.row(i).t())) * 
                (Sit.each_col() - zc.row(i).t()) * dl.subvec(i,t).t();
        }
        ht[t] = hhat;
    }

    return(ht);
}

En mi sistema, se trata de un factor de dos más rápido que su código. Sin embargo, es posible que sea aún más rápido.

2
Ralf Stubner 21 dic. 2019 a las 23:41