I have written an IIR filter class with an internal filter state and an element-wise function float filter(float x). Naturally, because of the internal filter state, filter() can only be called on the samples in the order in which they are arranged in the sequence.
If I now apply that function to an Eigen Array, like
input.unaryExpr(filter);
will the order of execution of unaryExpr() always be strictly in order of the values in the array, or may there out-of-order execution or even parallelization happen?
Would it be more safe to write the loop explicitly, to make sure the order is always as expected?
Right now it seems to be working correct, but I can't find any explicit documentation on its behavior.
A unary expression, like all Eigen expressions, is really just a fancy functor that behaves like an Eigen::Matrix or Array, that allows some introspection, and can be asked for the values for individual row and column indices. The order in which these values are requested is determined by the object to which the values are assigned; essentially, the left side of the assignment. We can demonstrate this in 2D:
Eigen::ArrayXXi in(2, 3);
in <<
0, 1, 2,
3, 4, 5;
// just the identity function with some debut output
auto expr = in.unaryExpr([](int x) {
std::cout << x << ' ';
return x;
});
// assignment 1, simple copy
Eigen::ArrayXXi out1 = expr;
std::cout << '\n';
// assignment 2, copy to row major array
using RowMajorArray = Eigen::Array<
int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
RowMajorArray out2 = expr;
std::cout << '\n';
// assignment 3, copy to transposed output array
Eigen::ArrayXXi out3(3, 2);
out3.transpose() = expr;
The first assignment prints 0 3 1 4 2 5, the other two assignments print 0 1 2 3 4 5. The reason is that out1 uses the default column-major order. So the most efficient order to fill it with values is column 0 row 0, column 0 row 1, column 1, row 0, column 1 row 1, … The other two output arrays have transposed memory order, therefore use transposed assignment order.
Note also how there is nothing stopping me from evaluating twice. Of course this here is not regular use, but it can crop up in otherwise benign code, like this template here:
template<class Derived>
typename Eigen::ArrayBase<Derived>::PlainObject
derivative(const Eigen::ArrayBase<Derived>& in)
{
Eigen::Index n = in.cols() - 1;
return in.rightCols(n) - in.leftCols(n);
}
If I call derivative(in.unaryExpr(…)), it will print 0 1 3 4 1 2 4 5, evaluating most entries twice.
Long story short: I don't think it is a good idea to make expressions stateful like you want, at least in general. Ideally, expressions should be idempotent. But if it's the shortest code, you know what you are doing, and you add some warnings around it so that the next programmer also understands the limitations, it should be okay.
In general, Eigen does not parallelize automatically except for matrix multiplications and similar, known complex operations. Everything else is assumed to be small and fast and mostly bound by memory bandwidth. I can't imagine this assumption changing except in a new major release Eigen4. And even then I would not assume it to be the default behavior for custom expressions. The cost of starting and stopping an OpenMP parallel region is just too high to do without good knowledge that it is justified. And by definition, Eigen has no clue about the computational complexity of your custom expression.
In general, the order of evaluation will also be front-to-back in the in-memory storage order of the output, at least for simple cases. That's what's fast to compute and straightforward to program.
Personally, I'd just use a plain old std::transform or a loop to make the order of operations explicit. Iterators to 1D arrays and vectors are fast, same with colwise() iterators in 2D.
Eigen::ArrayXf input(1234), output(1234);
std::transform(std::begin(input), std::end(input), std::begin(output),
[state](float x) mutable { return state(x); });