Subsetting Eigen vectors and matrices with a vector of indices

Question

I'm trying to port a working Armadillo function to Eigen and am having an issue with RcppEigen vector and matrix subsetting.

Here's my function:

//[[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);
}

It's throwing three errors all on the line indicated above:

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> >&)'

Here's my RcppArmadillo equivalent (working):

//[[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);
}

I'm unsure why Eigen cannot subset using these indices. My implementation seems to be consistent with the Eigen Sub-Matrices documentation.

Any idea why this is throwing an error?

p.s. I've been able to use this same function with RcppArmadillo using the .submat and .elem functions with a uvec indices vector generated by arma::find. There is apparently no equivalent to arma::find in Eigen.

UPDATE I've found documentation directly on this, and I think we can expect support for non-contiguous subviews of Eigen matrices in the (near) future:

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

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

Answer 1

I may be reading the Eigen documentation differently: I do not think you can 'pick' elements from a matrix or vector by injecting an integer vector. If it did as you do above with nz then the simpler below would compile. But it doesn't. Meaning your very clever and very highly aggregate 'update' expression does not work.

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

There is additional documentation (but that is also from Eigen 3.4.*) suggesting something closer to what you use with Armadillo but I have not tried this.

Answer 2

Apparently this is a much-requested feature in Eigen and is coming in 4.0 (yay!)

As Dirk pointed out, there is no way to do this without copying data to a new matrix, but my microbenchmarking suggests the below Eigen functions do so about 20% faster than Armadillo .submat() and .subvec().

Here is the function for subsetting an object of 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));
}

And here is how I used that function in my nnls function. This implementation also shows how to subset a vector:

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