Rcpp操作矩阵和向量集锦

Posted at — Jan 7, 2016

收集和记录Rcpp或者RcppArmadillo操作矩阵和向量。

1. Rcpp

可以使用逻辑下标（LogicalVector）对向量和列表取值。

2. RcppArmadillo

基本类型是mat、vec（colvec）和rowvec。

属性
- 对于矩阵，行数：m.n_rows;；列数：m.n_cols;；维度：m.size();或size(m);。对于向量，元素数：v.n_elem;。
特殊向量或矩阵
- 全是1ones<mat>(3, 4);/vec(10, fill::ones);/；全是0zeros<vec>(10);/mat(3, 4, fill::zeros);；全是某个数mat a(4, 5); a.fill(123.4);。
- 连续向量，规定长度linspace<vec>(0, 5, 6);；连续向量，规定间距regspace<vec>(0, 2, 9);。
取值
- 对于向量，连续取值：v.subvec(1stIdx, lastIdx);；非连续，可以考虑使用find()函数，比如：v.elem(find(v > 0));。
- 对于矩阵，连续取值：m.col(Idx);/m.row(Idx);/m.cols(Idx);/m.rows(Idx);/m.submat(1stRowIdx, lastRowIdx, 1stColIdx, lastColIdx);；非连续，m.submat(vecRowIdx, vecColIdx);。
Rcpp矩阵转换为RcppArmadillo矩阵，可以避免拷贝矩阵，以提升效率，mat(ptr_aux_mem, n_rows, n_cols, copy_aux_mem = true, strict = false);。同样道理，可以转化向量。例如：

arma::mat TransferMatArma(Rcpp::NumericMatrix x, Rcpp::NumericVector y) {
    mat tx(x.begin(), x.nrow(), x.ncol(), false);
    vec ty(y.begin(), y.size(), false);
    return tx;
}

Rcpp::NumericVector TransferMatRcpp(arma::mat x, arma::vec y) {
    NumericMatrix tx(x.n_rows, x.n_cols, x.begin());
    NumericVector ty(y.begin(), y.end());
    return ty;
    
// 不要使用as<IntegerVector>(wrap(y))，会有内存泄露。
}

使用.each_col()/.each_row()/.for_each()替代apply()

arma::mat TestMat(arma::mat M, double a) {

  M.for_each([a](mat::elem_type& val) {
      val = val > 0 ? val : a;
    });

  M.each_row([](rowvec& r) {
      r /= r.max();
    });

  return M;
}

使用sum(M, 0);和sum(M, 1);分别替代colSums(M)和rowSums(M)。
使用join_rows()和join_cols()分别替代cbind()和rbind()。注意事项：连续合并列或行会造成性能下降。好的方法是：提前建立向量，然后将列或行依次填入。

3. bigmemory

bigmemory包提供了四种数据类型的矩阵，即double（默认）、integer、short和char。对于big.matrix对象pMat，四种类型通过通过MatrixAccessor<double> macc(*pMat)、MatrixAccessor<int> macc(*pMat)、MatrixAccessor<short> macc(*pMat)和MatrixAccessor<char> macc(*pMat)提取元素。pMat有三种属性nrow()、ncol()和matrix_type()可以使用。一下代码示例展示了将big.matrix转换为matrix：

#include <Rcpp.h>
// [[Rcpp::depends(BH, bigmemory)]]
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>

using namespace Rcpp;


// [[Rcpp::export]]
Rcpp::NumericMatrix TestBigMat(XPtr<BigMatrix> pMat) {

  MatrixAccessor<int> macc(*pMat);

  int n = pMat->nrow();
  int m = pMat->ncol();

  NumericMatrix resMat(n, m);

  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      resMat(i, j) = macc[j][i];
    }
  }

  return resMat;
}

注意事项：

获取元素为列-行形式，因为big.matrix按照列存储矩阵。例如macc[j][i]表示i-1行的j-1元素。
调用函数使用big.matrix的地址，例如TestBigMat(bigmat@address)。

同样，RcppArmadillo也能与bigmemory结合，例如：

#include <RcppArmadillo.h>
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>

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

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
Rcpp::NumericMatrix TestBigMatArma(SEXP pMat) {

  XPtr<BigMatrix> xpMat(pMat);

  Mat<int> macc = Mat<int>((int *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false);

  int n = xpMat->nrow();
  int m = xpMat->ncol();

  NumericMatrix resMat(n, m);

  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      resMat(i, j) = macc(i, j);
    }
  }

  return resMat;
}

参考网址

更新记录

2021年12月08日

Yulong Niu

个人博客