About the use of .shed_row/.shed_col in RcppArmadillo

Question

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!!!

Answer 1

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.

Answer 2

The docs are good here. The rows / columns are removed in place. So after changing the return type and the type of the matrix multiplication and remembering the indices start at zero in c++, it is enough to do

#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)
*/

---

Re your update; another strategy is to select the columns of the matrices to retain, rather than removing a column. This is done by creating an sequence along the columns with arma::regspace and then retaining the columns which are not to be removed with 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))
*/