用索引向量对特征向量和矩阵进行子集化

Question

我正在尝试将一个可用的 Armadillo 函数移植到 Eigen，但我遇到了 RcppEigen 向量和矩阵子集设置的问题。

这是我的函数：

//[[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Eigen;

// [[Rcpp::export]]
Eigen::VectorXd fastnnls_eigen(const Eigen::MatrixXd& a, const Eigen::VectorXd& b, int maxit = 50) {
  Eigen::VectorXd x = (a).llt().solve((b));
  
  while(maxit-- > 0){

    // find values in x greater than zero
    // set values less than zero to zero
    bool x_is_nonneg = true;
    std::vector<int> nz;
    for(int i = 0; i < x.size(); ++i){
      if(x(i) > 0){
        nz.push_back(i);
      }
      else if(x(i) < 0) {
        x_is_nonneg = false; 
        x(i) = 0;
      }
    }
    if(x_is_nonneg) break;

    // update x with solutions from only indices given in "nz"
    x(nz) = a(nz, nz).llt().solve((b(nz)));   // *************ERROR ON THIS LINE
  }
  return(x);
}

它在上面指出的行中抛出了三个错误：

no match for call to '(Eigen::VectorXd {aka Eigen::Matrix<double, -1, 1>}) (std::vector<int, std::allocator<int> >&)'
no match for call to '(Eigen::MatrixXd {aka Eigen::Matrix<double, -1, 1>}) (std::vector<int, std::allocator<int> >&)'
no match for call to '(Eigen::VectorXd {aka Eigen::Matrix<double, -1, 1>}) (std::vector<int, std::allocator<int> >&)'

这是我的 RcppArmadillo 等价物（工作）：

//[[Rcpp::export]]
arma::vec fastnnls(const arma::mat& a, const arma::vec& b) {
  arma::vec x = arma::solve(a, b, arma::solve_opts::likely_sympd + arma::solve_opts::fast);
  while (arma::any(x < 0)) {
    arma::uvec nz = arma::find(x > 0);
    x.fill(0);
    x.elem(nz) = arma::solve(a.submat(nz, nz), b.elem(nz), arma::solve_opts::likely_sympd + arma::solve_opts::fast);
  }
  return(x);
}

我不确定为什么 Eigen 不能使用这些索引进行子集化。我的实现似乎与特征子矩阵一致 documentation.

知道为什么会抛出错误吗？

p.s。我已经能够使用 RcppArmadillo 使用 .submat 和 .elem 函数以及 arma::find 生成的 uvec 索引向量来使用相同的函数。在 Eigen 中显然没有 arma::find 等价物。

更新我直接找到了关于此的文档，我认为我们可以期待在（不久的）将来支持 Eigen 矩阵的非连续子视图：

https://gitlab.com/libeigen/eigen/-/issues/329

http://eigen.tuxfamily.org/dox-devel/TopicCustomizing_NullaryExpr.html#title1

Answer 1

我阅读 Eigen 文档的方式可能有所不同：我认为您无法通过注入整数向量来 'pick' 矩阵或向量中的元素。如果它像上面那样用 nz 做，那么下面的更简单的将编译。但事实并非如此。这意味着您非常聪明且高度聚合的 'update' 表达式不起作用。

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

// [[Rcpp::export]]
Eigen::VectorXd demoSubset(const Eigen::VectorXd& b, std::vector<int> p) {
  return b(p);
}

/*** R
demoSubset(as.numeric(1:10), c(2L,4L,8L))
*/

additional documentation（但这也是来自 Eigen 3.4.*）建议使用更接近于 Armadillo 的方法，但我还没有尝试过。

Answer 2

显然，这是 Eigen 中一个非常受欢迎的功能，即将在 4.0 中推出（耶！）

正如 Dirk 指出的那样，如果不将数据复制到新矩阵就无法做到这一点，但我的微基准测试表明下面的 Eigen 函数比 Armadillo .submat() 和 [=13 快 20% =].

这是对 class Eigen::MatrixBase 的对象进行子集化的函数：

http://eigen.tuxfamily.org/dox-devel/TopicCustomizing_NullaryExpr.html#title1

template<class ArgType, class RowIndexType, class ColIndexType>
class indexing_functor {
  const ArgType& m_arg;
  const RowIndexType& m_rowIndices;
  const ColIndexType& m_colIndices;
public:
  typedef Matrix<typename ArgType::Scalar,
    RowIndexType::SizeAtCompileTime,
    ColIndexType::SizeAtCompileTime,
    ArgType::Flags& RowMajorBit ? RowMajor : ColMajor,
    RowIndexType::MaxSizeAtCompileTime,
    ColIndexType::MaxSizeAtCompileTime> MatrixType;

  indexing_functor(const ArgType& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
    : m_arg(arg), m_rowIndices(row_indices), m_colIndices(col_indices) {}

  const typename ArgType::Scalar& operator() (Index row, Index col) const {
    return m_arg(m_rowIndices[row], m_colIndices[col]);
  }
};

template <class ArgType, class RowIndexType, class ColIndexType>
CwiseNullaryOp<indexing_functor<ArgType, RowIndexType, ColIndexType>, typename indexing_functor<ArgType, RowIndexType, ColIndexType>::MatrixType>
submat(const Eigen::MatrixBase<ArgType>& arg, const RowIndexType& row_indices, const ColIndexType& col_indices) {
  typedef indexing_functor<ArgType, RowIndexType, ColIndexType> Func;
  typedef typename Func::MatrixType MatrixType;
  return MatrixType::NullaryExpr(row_indices.size(), col_indices.size(), Func(arg.derived(), row_indices, col_indices));
}

下面是我在 nnls 函数中使用该函数的方式。此实现还展示了如何对向量进行子集化：

template<typename T, typename Derived>
Eigen::Matrix<T, -1, 1> nnls(const Eigen::Matrix<T, -1, -1>& a, const typename Eigen::MatrixBase<Derived>& b) {
  
  // initialize with unconstrained least squares solution
  Eigen::Matrix<T, -1, 1> x = a.llt().solve(b);
  while ((x.array() < 0).any()) {
    // get indices in "x" greater than zero (the "feasible set")
    Eigen::VectorXi gtz_ind = find_gtz(x);

    // subset "a" and "b" to those indices in the feasible set
    Eigen::Matrix<T, -1, 1> bsub(gtz_ind.size());
    for (unsigned int i = 0; i < gtz_ind.size(); ++i) bsub(i) = b(gtz_ind(i));
    Eigen::Matrix<T, -1, -1> asub = submat(a, gtz_ind, gtz_ind);

    // solve for those indices in "x"
    Eigen::Matrix<T, -1, 1> xsub = asub.llt().solve(bsub);
    x.setZero();
    for (unsigned int i = 0; i < nz.size(); ++i) x(nz(i)) = xsub(i);
  }
  return x;
}

用索引向量对特征向量和矩阵进行子集化

Subsetting Eigen vectors and matrices with a vector of indices

c++

rcpp

eigen

rcppeigen