关于RcppArmadillo中.shed_row/.shed_col的使用

Question

我现在正在尝试将 R 代码转换为 Rcpp 代码。
R代码是

hh <- function(A,B,j){

    aa <- A[j,-j] %*% B[j,-j] ## A and B are n*m matrixes
    return(aa)
}

> set.seed(123)
> A <- matrix(runif(30),5,6)
> B <- matrix(rnorm(30),5,6)
> j <- 2
> hh(A,B,j)
>           [,1]
> [1,] 0.9702774

我的 Rcpp 代码是

#include <RcppArmadillo.h>

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

double hh(arma::mat A, arma::mat B, arma::uword j){
  
  arma::mat Bj = B.shed_col(j);  /* error occurs */
  
  arma::mat Ak = A.row(j);
  
  double aa = Ak.shed_col(j) * arma::trans(Bj.row(j));  /* error occurs */
  
  return aa;
  
}

错误应该是关于.shed_row/.shed_col的使用。我用 Google 搜索了 .shed_row，但是，还没有解决我在这里遇到的问题的想法。你有什么主意吗？提前致谢！

进一步更新：
现在我们考虑在函数的for循环中使用.shed_row/.shed_col。
具体来说，我的Rcpp代码如下

#include <RcppArmadillo.h>

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

arma::mat ggcpp(arma::mat A, arma::mat B){
  
  /*we assume here A and B are n*n square matrices.*/
  
  int Ac = A.n_cols;
  int Bc = B.n_cols;
  arma::mat C(Ac,Bc);

  for(arma::uword i = 0; i < Ac; i++){
    
    A.shed_col(i);
    
    for(arma::uword j = 0; j < Bc; j++){
      
      B.shed_col(j);
      
      C(i,j) = arma::as_scalar(A.row(i) * B.row(j).t()); 
    }
  }
 
 return C; 
}

等效的R代码如下

gg <- function(A,B){
  
  ac <- ncol(A)
  bc <- ncol(B)
  C <- matrix(NA,ac,bc)
  
  for(i in 1:ac){
    
    for(j in 1:bc){
      
      C[i,j] <- A[i,-i] %*% B[j,-j]
    }
  }
  
  return(C)
}

我的 R 代码有效。它已经过测试。但是，我在使用 Rcpp 代码时遇到了麻烦。
试了很多方法，主要遇到两个错误：

ggcpp(a1, a2) 错误：as_scalar()：尺寸不兼容
ggcpp(a1, a2) 错误：Mat::shed_col()：索引越界

这里，a1和a2是两个随机生成的6*6矩阵

你有什么想法吗？赞赏！！！

Answer 1

docs这里不错。行/列被原地移除。因此，在更改 return 类型和矩阵乘法的类型并记住索引在 C++ 中从零开始后，就足够了

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat hhcpp(arma::mat A, arma::mat B, arma::uword j){ 

  A.shed_col(j); B.shed_col(j);
  arma::mat aa = A.row(j) * B.row(j).t(); 
  return aa;
}

/***R
 set.seed(123)
 A <- matrix(runif(30),5,6)
 B <- matrix(rnorm(30),5,6)
 j <- 2
 
 hh <- function(A,B,j){
    aa <- A[j,-j] %*% B[j,-j] ## A and B are n*m matrixes
    return(aa)
 }
 hh(A,B,j)
 
 hhcpp(A, B, j-1)
*/

回复你的更新；另一种策略是 select 保留矩阵的列，而不是删除列。这是通过使用 arma::regspace 沿列创建一个序列，然后保留不使用 find.

删除的列来完成的

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

arma::mat ggcpp(arma::mat A, arma::mat B){
  int Ac = A.n_cols;
  int Bc = B.n_cols;
  arma::mat C(Ac,Bc);
  
  arma::vec aRange = arma::regspace<arma::vec>(0, Ac-1);
  arma::vec bRange = arma::regspace<arma::vec>(0, Bc-1);

  for(int i = 0; i < Ac; i++){

    arma::mat Atemp = A.submat(i* arma::ones<arma::uvec>(1), find(aRange != i)) ;

    for(int j = 0; j < Bc; j++){

      arma::mat Btemp = B.submat(j* arma::ones<arma::uvec>(1), find(bRange != j)) ;

      C(i,j) = arma::as_scalar(Atemp * Btemp.t());
    }
   }
 
 return C; 
}

/***R
 # Square matrices only!!
 set.seed(123)
 A <- matrix(runif(30),5,5)
 B <- matrix(rnorm(30),5,5)
 j <- 2

 gg <- function(A,B){
  ac <- ncol(A)
  bc <- ncol(B)
  C <- matrix(NA,ac,bc)
  for(i in 1:ac){
    for(j in 1:bc){
      C[i,j] <- A[i,-i] %*% B[j,-j]
    }
  }
  return(C)
 }
 
 all.equal(gg(A,B), ggcpp(A, B))
*/

Answer 2

回答进一步更新：

#include <RcppArmadillo.h>

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

  Rcpp::NumericMatrix row_erase (Rcpp::NumericMatrix& x, Rcpp::IntegerVector& rowID) {
  rowID = rowID.sort();

  Rcpp::NumericMatrix x2(Rcpp::Dimension(x.nrow()- rowID.size(), x.ncol()));
  int iter = 0; 
  int del = 1; // to count deleted elements
  for (int i = 0; i < x.nrow(); i++) {
    if (i != rowID[del - 1]) {
      x2.row(iter) = x.row(i);
      iter++;
    } else {
      del++;
    }
  }
  return x2;
}


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

Rcpp::NumericMatrix col_erase (Rcpp::NumericMatrix& x, Rcpp::IntegerVector& colID) {
  colID = colID.sort();

  Rcpp::NumericMatrix x2(Rcpp::Dimension(x.nrow(), x.ncol()- colID.size()));
  int iter = 0; 
  int del = 1; 
  for (int i = 0; i < x.ncol(); i++) {
    if (i != colID[del - 1]) {
      x2.column(iter) = x.column(i);
      iter++;
    } else {
      del++;
    }
  }
  return x2;
}

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

arma::mat ggcpp(arma::mat A, arma::mat B){
    
    Rcpp::NumericMatrix AA = Rcpp::wrap(A);
    Rcpp::NumericMatrix BB = Rcpp::wrap(B);
    Rcpp::NumericMatrix AAi;
    Rcpp::NumericMatrix BBj;
    Rcpp::IntegerVector AiV;
    Rcpp::IntegerVector BjV;
    arma::mat Ai;
    arma::mat Bj;
    unsigned int Ac = A.n_cols;
    unsigned int Bc = B.n_cols;
    arma::mat C(Ac,Bc);
    
    for(arma::uword i = 0; i < Ac; i++){
        
        AiV = {i};
        AAi = col_erase(AA,AiV);
        Ai = Rcpp::as<arma::mat>(AAi);
        
        for(arma::uword j = 0; j < Bc; j++){
            
            BjV = {j};
            BBj = col_erase(BB,BjV);
            Bj = Rcpp::as<arma::mat>(BBj);
            C(i,j) = arma::as_scalar(Ai.row(i) * Bj.row(j).t());    
        }
    }
    
    return C;
}

注意：row_erase和col_erase是从借来的。

关于RcppArmadillo中.shed_row/.shed_col的使用

About the use of .shed_row/.shed_col in RcppArmadillo

r

rcpp

rcpparmadillo