使用C ++在R中对复杂矩阵进行乘法运算

Question

使用C ++在R中对复杂矩阵进行乘法运算

3

假设A是一个复矩阵。我想要高效地计算在R中的A%*%Conj(t(A))乘积。据我所知，使用C++可以显著加快速度，所以这就是我想要做的。

我有以下代码用于实矩阵，可以在R中使用。

library(Rcpp); 
library(inline); 
library(RcppEigen);

crossprodCpp <- '
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::Lower;

const Map<MatrixXd> A(as<Map<MatrixXd> >(AA));
const int   m(A.rows());
MatrixXd    AAt(MatrixXd(m, m).setZero().selfadjointView<Lower>().rankUpdate(A));
return wrap(AAt);
'

fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
A<-matrix(rnorm(100^2),100)
all.equal(fcprd(A),tcrossprod(A))

fcprd(A) 在我的笔记本电脑上比 tcrossprod(A) 运行得快得多。这是我在 A<-matrix(rnorm(1000^2),1000) 上得到的结果:

 microbenchmark::microbenchmark('tcrossprod(A)'=tcrossprod(A),'A%*%t(A)'=A%*%t(A),fcprd=fcprd(A))
 Unit: milliseconds
          expr       min       lq     mean   median       uq      max neval
 tcrossprod(A) 428.06452 435.9700 468.9323 448.8168 504.2628 618.7681   100
      A%*%t(A) 722.24053 736.6197 775.4814 767.7668 809.8356 903.8592   100
         fcprd  95.04678 100.0733 111.5021 103.6616 107.2551 197.4479   100

然而，这段代码仅适用于具有双精度条目的矩阵。我该如何修改此代码以使其适用于复杂矩阵？

我对编程了解非常有限，但我正在努力学习。任何帮助都将不胜感激！

- Cm7F7Bb

2个回答

1

这是一种将 Eigen::Map<Eigen::MatrixXcd> 对象绑定到 Rcpp 的方法。该解决方案适用于 R 包设置，但我不确定是否有一种简单的方法可以使用 inline 库组合它。

首先，您需要在 inst/include/mylib.h 中提供以下专业化，以便在 RcppExports.cpp 中包含此头文件：

#include <complex>

#include <Eigen/Core>
#include <Eigen/Dense>

#include <Rcpp.h>

namespace Rcpp   {
namespace traits {

  template<>
  class Exporter<Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> > > {
    using OUT = typename Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> >;
    const static int RTYPE = ::Rcpp::traits::r_sexptype_traits<std::complex<double>>::rtype;
    Rcpp::Vector<RTYPE> vec;
    int d_ncol, d_nrow;

  public:
    Exporter(SEXP x)
      : vec(x), d_ncol(1)
      , d_nrow(Rf_xlength(x)) {
      if (TYPEOF(x) != RTYPE)
        throw std::invalid_argument("Wrong R type for mapped matrix");
      if (::Rf_isMatrix(x)) {
        int* dims = INTEGER(::Rf_getAttrib(x, R_DimSymbol));
        d_nrow = dims[0];
        d_ncol = dims[1];
      }
    }
    OUT get() { return OUT(reinterpret_cast<std::complex<double>*>(vec.begin()), d_nrow, d_ncol); }
  };
}}

唯一的区别是在最后一行使用了 reinterpret_cast 的未专门化的 Exporter，而 std::complex 和 Rcomplex 都具有 C99 复合兼容类型，因此无论实现方式如何，它们都应该具有相同的内存布局。

总结一下，您现在可以创建以下函数：

// [[Rcpp::export]]
Eigen::MatrixXd selfadj_mult(const Eigen::Map<Eigen::MatrixXcd>& mat) {
  Eigen::MatrixXd result = (mat * mat.adjoint()).real();
  return result;
}

然后在R中调用函数：

library(mylib)
library(microbenchmark)

N <- 1000
A <- matrix(complex(real = rnorm(N * N), imaginary = rnorm(N * N)), N)

microbenchmark::microbenchmark(
    base  = A %*% Conj(t(A))
  , eigen = mylib::selfadj_mult(A)
  , times = 100L
)

代码在CentOS7 / GCC83上进行编译，使用-O3 -DNDEBUG -flto -march = generic。 R已从源代码构建，并使用完全相同的编译器/标志（使用默认的BLAS绑定）。结果如下：

Unit: seconds
  expr       min        lq      mean    median        uq       max neval
  base 2.9030320 2.9045865 2.9097162 2.9053835 2.9093232 2.9614318   100
 eigen 1.1978697 1.2004888 1.2134219 1.2031046 1.2057647 1.3035751   100

- Jean-Mathieu Vermosen

关于 library(inline)，我的理解是 cxxfunction 的 otherdefs 参数可以用来提供前向声明，因此我认为上面的特化也可以以同样的方式传递。 - Jean-Mathieu Vermosen

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ralf Stubner · Accepted Answer

Eigen库还支持通过使用复杂条目。因此，如果您将MatrixXd替换为MatrixXcd，则原则上应该可以工作。但是，这不会编译，可能是因为没有使用Map的复杂矩阵的as函数（参见https://github.com/RcppCore/RcppEigen/blob/master/inst/unitTests/runit.RcppEigen.R#L205）。需要as函数以在R数据类型和C++/Eigen数据类型之间进行转换（参见http://dirk.eddelbuettel.com/code/rcpp/Rcpp-extending.pdf）。如果您不使用Map，则可以使用以下内容：

crossprodCpp <- '
using Eigen::MatrixXcd;
using Eigen::Lower;

const MatrixXcd A(as<MatrixXcd>(AA));
const int   m(A.rows());
MatrixXcd    AAt(MatrixXcd(m, m).setZero().selfadjointView<Lower>().rankUpdate(A));
return wrap(AAt);
'

fcprd <- inline::cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen")
N <- 100
A <- matrix(complex(real = rnorm(N), imaginary = rnorm(N)), N)
all.equal(fcprd(A), A %*% Conj(t(A)))

然而，在我的测试中，这比基本的R版本要慢：

N <- 1000
A <- matrix(complex(real = rnorm(N * N), imaginary = rnorm(N * N)), N)
all.equal(fcprd(A), A %*% Conj(t(A)))
#> [1] TRUE
microbenchmark::microbenchmark(base = A %*% Conj(t(A)), eigen = fcprd(A))
#> Unit: milliseconds
#>   expr      min       lq     mean   median       uq      max neval
#>   base 111.6512 124.4490 145.7583 140.9199 160.3420 241.8986   100
#>  eigen 453.6702 501.5419 535.0192 537.2925 564.8746 628.4999   100

请注意，R中的矩阵乘法是通过BLAS完成的。然而，R使用的默认BLAS实现速度不太快。提高R性能的一种方法是使用优化的BLAS库，参见https://csgillespie.github.io/efficientR/set-up.html#blas-and-alternative-r-interpreters。

或者，如果您有完整的BLAS，则可以使用BLAS函数zherk。非常粗略：

dyn.load("/usr/lib/libblas.so")

zherk <- function(a, uplo = 'u', trans = 'n') {
    n <- nrow(a)
    k <- ncol(a)
    c <- matrix(complex(real = 0, imaginary = 0), nrow = n, ncol = n)
    z <- .Fortran("zherk",
             uplo = as.character(uplo),
             trans = as.character(trans),
             n = as.integer(n),
             k = as.integer(k),
             alpha = as.double(1),
             a = as.complex(a),
             lda = as.integer(n),
             beta = as.double(0),
             c = as.complex(c),
             ldc = as.integer(n))
    matrix(z$c, nrow = n, ncol = n)
}

N <- 2
A <- matrix(complex(real = rnorm(N * N), imaginary = rnorm(N * N)), N)
zherk(A, uplo = "l") - A %*% Conj(t(A))

请注意，这只填充了上（或下）三角部分，但速度非常快：

microbenchmark::microbenchmark(base = A %*% Conj(t(A)), blas = zherk(A))
#> Unit: milliseconds
#>  expr      min        lq      mean    median       uq      max neval
#>  base 112.5588 117.12531 146.10026 138.37565 167.6811 282.3564   100
#>  blas  66.9541  70.12438  91.44617  82.74522 108.4979 188.3728   100