Создать список RcppArmadillo матриц

Глубоко внутри алгоритма MCMC мне нужно умножить предоставленный пользователем список матриц с вектором, то есть следующий фрагмент Rcpp и RcppArmadillo кода вызывается несколько раз в MCMC итерации:Создать список RcppArmadillo матриц

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); 
}

Предоставленный пользователем список Mats остается фиксированным во время MCMC, вектор y изменений на каждой итерации. Эффективность имеет первостепенное значение, и я пытаюсь проверить, могу ли я ускорить код, не переустанавливая элементы Mats на arma :: mat, которые много раз (это нужно делать только один раз). Я попробовал следующий подход

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); 
}

, а затем

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); 
}

, но это не похоже на работу. Любые указатели альтернативных/лучших решений очень приветствуются.

источник

2016-08-08 user2692802

«Список» - это R (совместимый тип Rcpp), который может хранить другие R (совместимые Rcpp) типы (и те, для которых у нас есть конвертеры) - поэтому вам может потребоваться 'wrap()' your' arma: : mat' сначала, чтобы помочь компилятору с переменной, совместимой с SEXP. –

Спасибо за указатель. Итак, если я правильно понимаю, 'Rcpp :: List' преобразует матрицы Armadillo в простые Rcpp-матрицы. Есть ли способ создать «список» армадильо-матриц, поэтому я могу размножить их с помощью вектора Армадильо ** без ** необходимости конвертировать их обратно? – user2692802

Взгляните на документацию Armadillo; Классы полей могут быть тем, что вы хотите. Если вы хотите вернуть данные в R, вам нужны типы Rcpp, о которых знают классы R и Rcpp. Это включает типы Armadillo через RcppArmadillo. Изучите несколько примеров ... –

Хорошо, я в основном писал ответ в комментарии, но потом мне пришло в голову, что мы уже предоставляют рабочий пример в заглушке, созданной RcppArmadillo.package.skeleton():

// [[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); 
}

Это возвращает список внешний продукт (а матрица) и скалярное произведение данного вектора.

Что касается быстрого, а что нет: рекомендую не догадываться, а скорее профиль и измерять как можно больше. Моя склонность состояла бы в том, чтобы сделать больше (автономный) код на C++ в Armadillo и только вернуться с самого конца, сводя к минимуму конверсии.

источник

2016-08-08 20:31:56

ответ

Смежные вопросы