Search code examples
rrcpprcpparmadillo

About the use of .shed_row/.shed_col in RcppArmadillo


I am now trying to convert R code to Rcpp code.
The R code is

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

My Rcpp code is

#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;
  
}

The error should be about the use of .shed_row/.shed_col. I have Googled .shed_row, however, did not have an idea yet to address the issue I encountered here. Do you have any idea? Thank you in advance!

Further Update:
Now we consider using .shed_row/.shed_col in for-loop in the function.
Specifically, my Rcpp code is the following

#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; 
}

The equivalent R code is the following

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)
}

My R code works. It has been tested. However, I am having troubles with the Rcpp code.
I tried many ways and mainly encountered two errors :

  1. Error in ggcpp(a1, a2) : as_scalar(): incompatible dimensions
  2. Error in ggcpp(a1, a2) : Mat::shed_col(): index out of bounds

Here, a1 and a2 are two randomly generated 6*6 matrices.

Do you have any idea? Appreciated!!!


Solution

  • Answer to Further Update:

    #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;
    }
    
    

    Note: row_erase and col_erase are borrowed from here.