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
And here is how I used that function in my nnls function. This implementation also shows how to subset a vector: