Rcpp: removing a column and row from a matrix

Question

I am trying to create a function that takes a matrix, nxp and index e, and returns the submatrix obtained by excluding the e-th column and e-row from X. I thought that the easiest way was to create a matrix n -1 xp-1 and insert angles around the cross formed in it by the e-th row and column. The same approach works fine with EigenRcpp using the Corner syntax. Rcpp doesn't seem to like assigning to ranges. I get the following error message:

error: non-static reference member 'Rcpp::SubMatrix<14>::MATRIX& Rcpp::SubMatrix<14>::m', can't use default assignment operator

I will copy the code below.

My questions: How to assign values ​​to a matrix block? Is there a better way to do this?

 #include <Rcpp.h> using namespace Rcpp; NumericMatrix elimMat(NumericMatrix X, int e){ int p = X.cols() - 1; int n = X.rows() - 1; int f = e - 1; NumericMatrix M = NumericMatrix(n, p); M(Range(0, f - 1), Range(0, f - 1)) = X(Range(0, f - 1), Range(0, f - 1)); //TopLeft same M(Range(0, f - 1), Range(f, p - 1)) = X(Range(0, f - 1), Range(f + 1, p)); //TopRight M(Range(f, n - 1), Range(0, f - 1)) = X(Range(f + 1, n), Range(0, f - 1)); //BottomLeft M(Range(f, n - 1), Range(f, p - 1)) = X(Range(f + 1, n), Range(f + 1, p)); //BottomRight return M; } 

Thanks for any help, Marco

Answer 1

Regarding this error

error: non-static reference element 'Rcpp :: SubMatrix <14> :: MATRIX & Rcpp :: SubMatrix <14> :: m, cannot use the default assignment operator

I do not think that you can (or anyone on the client side) do anything, since the problem comes from the definition of the SubMatrix class. According to the compiler,

 private: MATRIX& m ; vec_iterator iter ; int m_nr, nc, nr ; 

m must be a member of static , and it is not, therefore, an error.

---

Although I cannot give you a general approach for using Range assignments in this way, here is one possible approach to your specific question about “deleting” the NumericMatrix section:

 #include <Rcpp.h> // [[Rcpp::export]] Rcpp::NumericMatrix crossout(Rcpp::NumericMatrix X, int e) { Rcpp::NumericMatrix result = Rcpp::no_init_matrix(X.nrow() - 1, X.ncol() - 1); Rcpp::NumericMatrix::iterator src = X.begin(), dst = result.begin(); while (src != X.end()) { if (((src - X.begin()) % X.nrow()) != e && ((src - X.begin()) / X.nrow()) != e) { *dst = *src; ++dst; } ++src; } return result; } 

---

 /*** R M <- matrix(1:9 + .5, nrow = 3) R> all.equal(M[c(1, 3), c(1, 3)], crossout(M, 1)) #[1] TRUE R> crossout(M, 1) # [,1] [,2] #[1,] 1.5 7.5 #[2,] 3.5 9.5 M2 <- matrix(1:20 + .5, nrow = 4) R> all.equal(M2[c(1:2, 4), c(1:2, 4:5)], crossout(M2, 2)) #[1] TRUE R> crossout(M2, 2) # [,1] [,2] [,3] [,4] #[1,] 1.5 5.5 13.5 17.5 #[2,] 2.5 6.5 14.5 18.5 #[3,] 4.5 8.5 16.5 20.5 R> any(unique(c(M2[3,], M2[,3])) %in% crossout(M2, 2)) #[1] FALSE */ 

---

Above, I copy the current value of the "source" iterator to the destination iterator if the source position of the iterator does not fall into the column or row specified by e . The row index is obtained as (src - X.begin()) % X.nrow() , and the column index is obtained as (src - X.begin()) / X.nrow() .

---

Edit: As for your comment below, it would be pretty simple to place non-square matrices / different indexes for row and column sizes:

 #include <Rcpp.h> // [[Rcpp::export]] Rcpp::NumericMatrix crossout2(Rcpp::NumericMatrix X, int r, int c) { if (r >= X.nrow() || c >= X.ncol()) { Rf_error("Invalid row or column index.\n"); } Rcpp::NumericMatrix result = Rcpp::no_init_matrix(X.nrow() - 1, X.ncol() - 1); Rcpp::NumericMatrix::iterator src = X.begin(), dst = result.begin(); while (src != X.end()) { if (((src - X.begin()) % X.nrow()) != r && ((src - X.begin()) / X.nrow()) != c) { *dst = *src; ++dst; } ++src; } return result; } 

---