More efficient way to get indices of a binary mask in Eigen3?

Question

I currently have a bool mask vector generated in Eigen. I would like to use this binary mask similar as in Python numpy, where depending on the True value, i get a sub-matrix or a sub-vector, where i can further do some calculations on these.

To achieve this in Eigen, i currently "convert" the mask vector into another vector containing the indices by simply iterating over the mask:

Eigen::Array<bool, Eigen::Dynamic, 1> mask = ... // E.G.: [0, 1, 1, 1, 0, 1];

Eigen::Array<uint32_t, Eigen::Dynamic, 1> mask_idcs(mask.count(), 1);
int z_idx = 0;
for (int z = 0; z < mask.rows(); z++) {
    if (mask(z)) {
        mask_idcs(z_idx++) = z;
     }
}
// do further calculations on vector(mask_idcs)
// E.G.:    vector(mask_idcs)*3 + another_vector

However, i want to further optimize this and am wondering if Eigen3 provides a more elegant solution for this, something like vector(from_bin_mask(mask)), which may benefit from the libraries optimization.

There are already some questions here in SO, but none seems to answer this simple use-case (1, 2). Some refer to the select-function, which returns an equally sized vector/matrix/array, but i want to discard elements via a mask and only work further with a smaller vector/matrix/array.

Is there a way to do this in a more elegant way? Can this be optimized otherwise?

(I am using the Eigen::Array-type since most of the calculations are element-wise in my use-case)

Answer 1

As far as I'm aware, there is no "out of the shelf" solution using Eigen's methods. However it is interesting to notice that (at least for Eigen versions greater or equal than 3.4.0), you can using a std::vector<int> for indexing (see this section). Therefore the code you've written could simplified to

Eigen::Array<bool, Eigen::Dynamic, 1> mask = ... // E.G.: [0, 1, 1, 1, 0, 1];

std::vector<int> mask_idcs;
for (int z = 0; z < mask.rows(); z++) {
    if (mask(z)) {
        mask_idcs.push_back(z);
    }
}
// do further calculations on vector(mask_idcs)
// E.G.:    vector(mask_idcs)*3 + another_vector

If you're using c++20, you could use an alternative implementation using std::ranges without using raw for-loops:

int const N = mask.size();
auto c = iota(0, N) | filter([&mask](auto const& i) { return mask[i]; });
auto masked_indices = std::vector(begin(c), end(c));
// ... Use it as vector(masked_indices) ...

I've implemented some minimal examples in compiler explorer in case you'd like to check out. I honestly wished there was a simpler way to initialize the std::vector from the raw range, but it's currently not so simple. Therefore I'd suggest you to wrap the code into a helper function, for example

auto filtered_indices(auto const& mask) // or as you've suggested from_bin_mask(auto const& mask)
{
    using std::ranges::begin;
    using std::ranges::end;
    using std::views::filter;
    using std::views::iota;

    int const N = mask.size();
    auto c = iota(0, N) | filter([&mask](auto const& i) { return mask[i]; });
    return std::vector(begin(c), end(c));
}

and then use it as, for example,

Eigen::ArrayXd F(5);
F << 0.0, 1.1548, 0.0, 0.0, 2.333;
auto mask = (F > 1e-15).eval();
auto D = (F(filtered_indices(mask)) + 3).eval();

It's not as clean as in numpy, but it's something :)

Answer 2

I have found another way, which seems to be more elegant then comparing each element if it equals to 0:

Eigen::SparseMatrix<bool> mask_sparse = mask.matrix().sparseView();
for (uint32_t k = 0; k<mask.outerSize(); ++k) {
    for (Eigen::SparseMatrix<bool>::InnerIterator it(mask_sparse, k); it; ++it) {
        std::cout << it.row() << std::endl;   // row index
        std::cout << it.col() << std::endl;   // col index
        // Do Stuff or built up an array
        }
    }

Here we can at least build up a vector (or multiple vectors, if we have more dimensions) and then later use it to "mask" a vector or matrix. (This is taken from the documentation).

So applied to this specific usecase, we simply do:

Eigen::Array<uint32_t, Eigen::Dynamic, 1> mask_idcs(mask.count(), 1);
Eigen::SparseVector<bool> mask_sparse = mask.matrix().sparseView();
int z_idx = 0;
for (Eigen::SparseVector<bool>::InnerIterator it(mask_sparse); it; ++it) {
    mask_idcs(z_idx++) = it.index()
}

// do Stuff like vector(mask_idcs)*3 + another_vector

However, i do not know which version is faster for large masks containing thousands of elements.