Rcpp:从矩阵中消除一列和一行 [英] Rcpp: Eliminating a column and a row from a matrix

问题描述

我正在尝试创建一个函数，该函数接受一个矩阵 nxp 和一个索引 e，并返回通过消除 获得的子矩阵第 e-th 列和来自 X 的 e-th 行.我认为最简单的方法是创建一个 n -1 xp - 1 矩阵并插入十字形周围的角由其中的e-th行和列形成.使用 Corner 语法的 EigenRcpp 也可以使用相同的方法.似乎 Rcpp 不喜欢给 Ranges 赋值.我收到以下错误消息:

I am trying to create a function that takes a matrix, n x p, and an index, e, and returns the sub-matrix obtained by eliminating the e-th column and the e-th row from X. I thought that the easiest thing was to create an n -1 x p - 1 matrix and insert the corners around the cross formed by the e-th row and column in it. The same approach works fine with EigenRcpp using the Corner syntax. It seems that Rcpp doesn't like assignment to Ranges. I get the following Error message:

错误:非静态引用成员 'Rcpp::SubMatrix<14>::MATRIX&Rcpp::SubMatrix<14>::m'，不能使用默认赋值运算符

我复制下面的代码.

我的问题:如何为矩阵块赋值?有没有更好的方法来做到这一点?

My questions: How do I assign values to a block of a matrix? 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;
}

感谢您的帮助，马可

推荐答案

关于这个错误，

错误:非静态引用成员‘Rcpp::SubMatrix<14>::MATRIX&Rcpp::SubMatrix<14>::m'，不能使用默认赋值运算符

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

我认为您(或客户端的任何人)对此无能为力，因为问题来自 SubMatrix 的类定义.根据编译器，

I don't think there is anything you (or anyone on the client side) can do about this since the problem is coming from the class definition of SubMatrix. According to the compiler,

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

m 应该是 static 成员，但不是，因此出现错误.

m should be a static member, and it's not, hence the error.

虽然我无法为您提供以这种方式使用 Range 分配的一般解决方法，但这里有一种可能的方法来解决您关于划掉"一段 的特定问题>数值矩阵:

Although I can't give you a general work-around for using Range assignments in this manner, here is one possible approach to your specific question about "crossing out" a section of a NumericMatrix:

#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

*/

<小时>

以上，如果源迭代器的位置不在 e 指定的列或行内，我会将源"迭代器的当前值复制到目标"迭代器.行索引获取为(src - X.begin()) % X.nrow()，列索引获取为(src - X.begin())/X.nrow().

---

Above, I'm copying the current value of the "source" iterator to the "destination" iterator if the source iterator's position does not fall within 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().

关于您在下面的评论，容纳行和列维度的非方阵/不同索引将非常简单:

Regarding your comment below, it would be pretty straightforward to accomodate non-square matrices / differing indices for the row and column dimensions:

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

这篇关于Rcpp:从矩阵中消除一列和一行的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！