2

In Rcpp/RcppArmadillo I want to do the following: From an n x n matrix A, I would like to extract a submatrix A[-j, -j] where j is a vector of indices: In R it can go like

A = matrix(1:16, 4, 4)
j = c(2, 3)
A[-j, -j]

Seems that this functionality is not available in Rcpp or RcppArmadillo - sorry if I have overlooked something. One approach in R is

pos = setdiff(1:nrow(A), j)
A[pos, pos]

That will carry over to RcppArmadillo, but it seems akward having to create the vector pos as the complement of j - and I am not sure how to do it efficiently.

Does anyone have an idea for an efficient implementation / or a piece of code to share?

user20650
  • 24,654
  • 5
  • 56
  • 91
Søren Højsgaard
  • 73
  • 3
  • The overarching problem is that `A[-i, -j]` is an R expression, governed by the R interpreter. Once you are in a line of code in `src/*` you are governed by the C (or C++, close enough) compiler where this means something else. There is not a lot we can do. You need a mediating layer and the provided answer is a pretty good one. – Dirk Eddelbuettel Aug 29 '22 at 11:50

1 Answers1

4

The armadillo docs have the function .shed which takes an argument(s) that "... contains the indices of rows/columns/slices to remove". From my reading to remove both rows and columns will take two calls of .shed().

Using your example

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

// [[Rcpp::export]]
arma::mat fun(arma::mat X, arma::uvec row, arma::uvec col) {
  X.shed_cols(col); // remove columns
  X.shed_rows(row); // remove rows
  return(X);
}


/***R
A = matrix(1:16, 4, 4)
j = c(2, 3)
A[-j, -j]

# minus one for zero indexing
fun(A, j-1, j-1)
*/
user20650
  • 24,654
  • 5
  • 56
  • 91