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

编程问答 2022-06-09

问题描述

我现在正在尝试将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),6)
> j <- 2
> hh(A,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代码时遇到了麻烦。
我尝试了很多方法，主要遇到两个错误：

  1. ggcpp（a1，a2）中的错误：as_scalar（）：不兼容的尺寸
  2. ggcpp（a1，a2）中的错误：Mat :: shed_col（）：索引超出范围

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

你有什么主意吗？感激！！！

解决方法

docs在这里很好。行/列被删除。因此，在更改返回类型和矩阵乘法的类型并记住索引在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),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,j)
 
 hhcpp(A,j-1)
*/

重新更新；另一种策略是选择要保留的矩阵列，而不是删除列。这是通过使用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)
 B <- matrix(rnorm(30),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))
*/
,

回答进一步更新：

#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是从here借来的。

rrrcpprcpparmadillo