Преобразование матрицы броненосца в собственную матрицуXd и наоборот

Как я могу преобразовать матрицу Armadillo в Eigen MatrixXd и наоборот?

У меня есть nu как arma::vec размера N, z как arma::mat размера N x 3. Я хочу вычислить матрицу P, такую ​​как запись P_ij

Pij=exp(nu(i) + nu(j) + z.row(j)*z.row(j)))

Таким образом, я использовал этот код

int N=z.n_rows;
mat P= exp(nu*ones(1,N) + one(N,1)*(nu.t()) + z*(z.t()));

Но вычисление занимает слишком много времени. В частности, для N = 50,000 время выполнения слишком велико.

Кажется, что использование Eigen может быть быстрее. Но моя матрица - Броненосец. Как я могу использовать собственные операции? Или как мне сделать эту операцию быстрее.


person Ari.stat    schedule 12.10.2017    source источник
comment
Чтобы ускорить умножение матриц в Armadillo, используйте OpenBLAS вместо стандартного BLAS. Это имеет огромное значение. Кроме того, ваш код использует функцию exp(), что может занять много времени. Вы можете ускорить его, включив OpenMP в Armadillo.   -  person hbrerkere    schedule 12.10.2017
comment
Танки Я тут не очень новичок. Пожалуйста, не могли бы вы привести простой пример?   -  person Ari.stat    schedule 12.10.2017
comment
Прочтите файл README.txt, прилагаемый к архиву Armadillo (страница скачать). См. также страницу FAQ Armadillo.   -  person hbrerkere    schedule 12.10.2017


Ответы (2)


Используя броненосца .memptr(), мы можем извлечь указатель памяти. Отсюда мы можем использовать Eigen Map<T>() для создания матрицы Eigen.

Теперь мы можем перейти от матрицы Eigen, используя .data() для извлечения точки в Структура памяти Эйгена. Затем, используя расширенные параметры конструктора arma::mat, мы можем создать броненосец матрица .

Например:

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

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

// [[Rcpp::export]]
Eigen::MatrixXd example_cast_eigen(arma::mat arma_A) {

  Eigen::MatrixXd eigen_B = Eigen::Map<Eigen::MatrixXd>(arma_A.memptr(),
                                                        arma_A.n_rows,
                                                        arma_A.n_cols);

  return eigen_B;
}

// [[Rcpp::export]]
arma::mat example_cast_arma(Eigen::MatrixXd eigen_A) {

  arma::mat arma_B = arma::mat(eigen_A.data(), eigen_A.rows(), eigen_A.cols(),
                               false, false);

  return arma_B;
}

/***R
(x = matrix(1:4, ncol = 2))
example_cast_eigen(x)
example_cast_arma(x)
*/

Полученные результаты:

(x = matrix(1:4, ncol = 2))
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

example_cast_eigen(x)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

example_cast_arma(x)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

Одно быстрое замечание: если вы используете функцию отображения Eigen, вы должны автоматически получить изменение в матрице Armadillo (и наоборот), например.

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

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

// [[Rcpp::export]]
void map_update(Eigen::MatrixXd eigen_A) {

  Rcpp::Rcout << "Eigen Matrix on Entry: " << std::endl << eigen_A << std::endl;

  arma::mat arma_B = arma::mat(eigen_A.data(), eigen_A.rows(), eigen_A.cols(),
                               false, false);

  arma_B(0, 0) = 10;
  arma_B(1, 1) = 20;

  Rcpp::Rcout << "Armadill Matrix after modification: " << std::endl << arma_B << std::endl;

  Rcpp::Rcout << "Eigen Matrix after modification: " << std::endl << eigen_A << std::endl;
}

Бежать:

map_update(x)

Выход:

Eigen Matrix on Entry: 
1 3
2 4

Armadill Matrix after modification: 
   10.0000    3.0000
    2.0000   20.0000

Eigen Matrix after modification: 
10  3
 2 20
person coatless    schedule 12.10.2017
comment
Ницца. Спасибо А как мне сконвертировать Эйгена в Броненосца - person Ari.stat; 12.10.2017
comment
Неплохой вопрос и хороший ответ. Может быть, дополнение к галерее? - person Dirk Eddelbuettel; 12.10.2017
comment
Да. Вероятно, в эти выходные с 3-мя другими сообщениями, которые мне нужно перенести... - person coatless; 12.10.2017
comment
@42 На вашем месте я бы использовал список рассылки rcpp-devel - person Dirk Eddelbuettel; 15.10.2017
comment
ХОРОШО. Я подпишусь. Кстати, что за Галерея? (Или, возможно, это уже есть в listinfo для rcpp-devel, и я увижу его через пару минут?) - person IRTFM; 15.10.2017
comment
Это широко цитируемый сайт здесь с исходным кодом здесь, в GH. - person Dirk Eddelbuettel; 15.10.2017
comment
Я признаю, что мне, вероятно, следовало погуглить, а не комментировать вопросы. Первый удар по очевидной стратегии. Итак, в покаянии я куплю вашу книгу. - person IRTFM; 15.10.2017
comment
@coatless Спасибо за ваш ответ. Ваше замечание очень важно, если я намерен обновить преобразованную матрицу. Я обнаружил, что если я не хочу, чтобы обновление затрагивало первую матрицу, я могу использовать arma(ptr_aux_mem, n_rows, n_cols, true, false) вместо arma(ptr_aux_mem, n_rows, n_cols, false, false). Но это для конвертации из Eigen в Armadillo. Знаете ли вы способ конвертировать из Armadillo в Eigen и разрешить обновление без изменения матрицы Armadillo? - person Ari.stat; 29.03.2018
comment
@coatless У меня проблема с этим. Первый элемент -т.е. at (0,0) - преобразования из Eigen в Armadillo не конвертируются правильно. Я не уверен, должен ли я публиковать новый вопрос. - person mkln; 01.10.2019
comment
@mkln Просто перезапустил выше, и все заработало. Дважды проверьте настройки компиляции, пакет и версии R. - person coatless; 01.10.2019
comment
pastebin.com/7AgTHPNf вот что я получаю. от R я также заставить его работать правильно - person mkln; 01.10.2019
comment
Ваш код создает новую копию при вызове вместо использования ссылки при вызове. Нужно сделать справочное обновление. Это отдельный вопрос. - person coatless; 01.10.2019
comment
То есть нужно использовать arma::mat matrixxd_to_armamat(Eigen::MatrixXd& eigen_A). - person coatless; 01.10.2019
comment
@coatless Я отправлю вопрос, я получу тот же результат, даже если воспользуюсь ссылкой - person mkln; 01.10.2019
comment
@coatless здесь: stackoverflow.com/questions/58176022/ - person mkln; 01.10.2019

Я просто трачу пару часов, пытаясь преобразовать разреженную матрицу Eigen в разреженную матрицу Armadillo, и я опубликую код здесь, если кто-то еще найдет необходимость сделать то же самое.

Я делал это, потому что не мог найти собственный решатель для разреженных комплексных матриц, а Armadillo был единственным, у кого он был, но остальная часть моего кода уже была написана в Eigen, поэтому мне пришлось выполнить преобразование.

#include <Eigen/Sparse>
#include <armadillo>
using namespace std;
using namespace arma;
int main() {
    auto matrixA = new SparseMatrix<complex<double>>(numCols, numRows); //your Eigen matrix

    /*
    SOME CODE TO FILL THE Eeigen MATRIX

    */
    //  now create a separate vectors for row indeces, first non-zero column element indeces and non-zero values
    //  why long long unsigned int, because armadilo will expect that type when constructing sparse matrix
    vector<long long unsigned int> rowind_vect((*matrixA).innerIndexPtr(),
                                               (*matrixA).innerIndexPtr() + (*matrixA).nonZeros());
    vector<long long unsigned int> colptr_vect((*matrixA).outerIndexPtr(),
                                               (*matrixA).outerIndexPtr() + (*matrixA).outerSize() + 1);
    vector<complex<double>> values_vect((*matrixA).valuePtr(),
                                        (*matrixA).valuePtr() + (*matrixA).nonZeros());

    //  you can delete the original matrixA to free up space
    delete matrixA;

    //new Armadillo vectors from std::vector, we set the flag copy_aux_mem=false, so we don't copy the values again
    cx_dvec values(values_vect.data(), values_vect.size(), false);
    uvec rowind(rowind_vect.data(), rowind_vect.size(), false);
    uvec colptr(colptr_vect.data(), colptr_vect.size(), false);

    //  now create Armadillo matrix from these vectors
    sp_cx_dmat arma_hamiltonian(rowind, colptr, values, numCols, numRows);

    //  you can delete the vectors here if you like to free up the space

    return 0; 

}    
person Aleksandar Bukva    schedule 06.04.2018