Stack Exchange
Subsetting Eigen vectors and matrices with a vector of indices
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

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


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.

Thanks Dirk!  I guess it's not possible then.  I'll explore memory casting between arma and eigen matrices (I see @coatless has a great example on this). The raw speed of Eigen is undeniable, but the elegance of Armadillo wins.

Yes. Both have had zero-copy constructors for a loooong time.  So you could subset with Arma, and then solve with Eigen.  Sounds complex.  Maybe a simpler 'reducer' helper to subset-by-copying-elements is easier?  Your call...

Just did some microbenchmarking. Memory casting takes nanoseconds or is indiscernible. This is much preferable to subsetting by copying elements, which would certainly take a bit longer. Armadillo makes castings easy, Eigen doesn't.

There is an element you may overlook: Armadillo likely copies to create a new subset matrix.  (There is likely no other way: vectors and matrices are contiguous chunks of memory which is why a) zero-copy ctor use is easy and fast b) likely why Eigen has no subsetter.) So I think you may not have compared fairly.  But I have no bone in this fight.  You had an error, we fixed it, if that is how you want to proceed and are merry with it, go have at it and enjoy!
