在 Rcpp Armadillo 中将 S4 对象转换为矩阵类型

Question

在我 RcppArmadillo 的项目中，我有一些矩阵（例如 mat A,B,C;）和一个 S4 对象，例如D（来自 R 中的外部函数）。由于我需要在这些矩阵和 D 之间进行一些计算，因此我想将 "D" 转换为 RcppArmadillo 中合适的数据类型，例如 arma::mat D 。可能吗？最好的方法是什么？是类似的代码：

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

// [[Rcpp::export]]
Rcpp::List func1(arma::mat A, arma::mat B){

  // Incoming
  Rcpp::List outcome;
  arma::mat rvecs;
  arma::vec rvals;
   Rcpp::Environment Matrix("package:Matrix"); 
   Rcpp::Function nearPD = Matrix["nearPD"];

 // Computation   
 Rcpp::List PD=nearPD(B);
 Rcpp::S4 D = PD["mat"];
 eig_sym(rvals, rvecs, D);
 arma::mat RI12_hat = rvecs * arma::diagmat(1.0/sqrt(rvals)) * rvecs.t();
  arma::mat diff = A - D;

  // Release results
  outcome = Rcpp::List::create(Rcpp::Named("rvals")    = rvals, 
                               Rcpp::Named("RI12_hat") = RI12_hat, 
                               Rcpp::Named("rvecs")    = rvecs);
  return outcome;   
}

其中D是classdpoMatrix的矩阵，计算出的正定矩阵，错误信息是“no matching for operator”。

Answer 1

不幸的是，dpoMatrix object yet. However, using the S4 class, we can extract out the necessary components and reuse the memory using armadillo's advanced ctor for matrices. First, we must understand the dpoMatrix object by either looking at its documentation for the underlying Matrix class 或构建示例没有 as<>() 或 wrap()。

考虑以下几点：

B <- matrix(1, 3,3); B[1,3] <- B[3,1] <- 0
n.B <- nearPD(B, corr=TRUE, do2eigen=FALSE)
str(n.B)

这给出：

List of 7
 $ mat        :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
  .. ..@ x       : num [1:9] 1 0.761 0.157 0.761 1 ...
  .. ..@ Dim     : int [1:2] 3 3
  .. ..@ Dimnames:List of 2
  .. .. ..$ : NULL
  .. .. ..$ : NULL
  .. ..@ uplo    : chr "U"
  .. ..@ factors : list()
 $ eigenvalues: num [1:3] 2.157 0.843 -0.679
 $ corr       : logi TRUE
 $ normF      : num 0.528
 $ iterations : num 18
 $ rel.tol    : num 6.48e-08
 $ converged  : logi TRUE
 - attr(*, "class")= chr "nearPD"

因此，我们可以从 x 槽中检索矩阵的值，并使用 .slot("name_here") 从 Dim 槽中检索矩阵的维度成员函数。

实施

话虽如此，我们将快速实现如下：

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

// [[Rcpp::export]]
Rcpp::List magic_func(arma::mat A, arma::mat B){

  // Incoming
  Rcpp::List outcome;
  arma::mat rvecs;
  arma::vec rvals;
  Rcpp::Environment Matrix("package:Matrix"); // Load the Matrix package in R!
  Rcpp::Function nearPD = Matrix["nearPD"];   // Extract nearPD() R function

  // Compute with R function an S4 object
  Rcpp::List PD = nearPD(B);
  Rcpp::S4 D_s4 = PD["mat"];

  // Convert the S4 object to an Armadillo matrix
  Rcpp::NumericVector temp = Rcpp::NumericVector(D_s4.slot("x"));
  Rcpp::NumericVector dims = D_s4.slot("Dim");

  // Advanced armadillo matrix ctor that reuses memory
  arma::mat D(temp.begin(), // pointer to NumericVector
              dims[0],      // Number of Rows
              dims[1],      // Number of Columns
              false,        // Avoid copying by disabling `copy_aux_mem`
              true);        // Bind memory by enabling `strict`

  // Computation
  eig_sym(rvals, rvecs, D);
  arma::mat RI12_hat = rvecs * arma::diagmat(1.0/sqrt(rvals)) * rvecs.t();
  arma::mat diff = A - D;

  // Return result
  outcome = Rcpp::List::create(Rcpp::Named("rvals")    = rvals, 
                               Rcpp::Named("RI12_hat") = RI12_hat, 
                               Rcpp::Named("rvecs")    = rvecs);
  return outcome;   
}

测试

代码：

set.seed(27)
A = matrix(round(rnorm(9),2), 3, 3)
A = A + t(A)

B = matrix(1, 3, 3); B[1,3] <- B[3,1] <- 0

magic_func(A, B)

结果：

$rvals
             [,1]
[1,] 2.414214e-08
[2,] 1.000000e+00
[3,] 2.414214e+00

$RI12_hat
          [,1]      [,2]      [,3]
[1,]  1609.647 -2275.222  1608.647
[2,] -2275.222  3218.293 -2275.222
[3,]  1608.647 -2275.222  1609.647

$rvecs
           [,1]          [,2]      [,3]
[1,] -0.5000000 -7.071068e-01 0.5000000
[2,]  0.7071068 -7.077672e-16 0.7071068
[3,] -0.5000000  7.071068e-01 0.5000000

在 Rcpp Armadillo 中将 S4 对象转换为矩阵类型

converting a S4 object to a matrix type in Rcpp Armadillo

r

armadillo

rcpp

实施

测试