Hello,

In my project, I have a big VectorXd of size n. I need to move a window of size m on it and to store all the different positions in a std::vector. It is easy if I create a std::vector<VectorXd>, but I don’t want to do that since I’ll fill the memory with a lot of useless copies. It’ll require approximately (n-m+1)*m*sizeof(double) bytes !

So, I created a class View deriving from Eigen::Map, and I added a template constructor that allows me to construct a View directly from any expression. Now I can easily put my Views in a std::vector. I’ve benchmarked a noticeable performance loss (certainly because of the dynamic strides), but it saves me a lot of memory and works very well with block expressions. Moreover, it also allows me to avoid template functions everywhere when I want to pass an Eigen type to a function.

My constructor needs the following methods to be available on any expression, to properly initialize the Map object: rowStride(), colStride() and data().

The problem is that they are not always accessible. For example, if my constructor gets a reversed vector:

Code: Select all

VectorXd vec(10);
MatrixXd mat(10,10);
double * a = vec.data(); // OK
double * b = vec.reverse().data();   // I was expected to get a pointer to the last element of vec, but the method doesn’t exist

I found a fix for this, by taking the address of the first element, like this:

Code: Select all

double * c = &vec.reverse()(0);   // ok it works, but in some cases the operator() doesn’t return a reference (but I forgot when… :s I think it was related to constness)

The same also applies for the expression returned by mat.diagonal()

But this is not the only problem I have:

Code: Select all

int d = vec.rowStride();   // OK
int f = mat.diagonal().rowStride(); // OK
int e = vec.reverse().rowStride();   // I was expected to get a negative stride, but this method returns void and is protected

To conclude, I have a last question: why isn’t it allowed to initialize an Eigen::Stride<Dynamic,Dynamic> with negative strides ? I think it could be useful to map a reversed matrix for example. (I imagine that it is not possible using static strides because of the ambiguity with the value of the keyword "Dynamic", but why should it also be forbidden with the dynamic strides?)

Thanks in advance for your help,

Fred

I think the main reason why we don't allow negative strides is that the benefit is relatively little (it mostly only benefits Reverse) for something that's a bit awkward since, as you mention, in the fixed-size case there is ambiguity with Dynamic=-1 and in general negative strides would break some naive code (when writing optimized low-level code it's quite easy to make the positive stride assumption without even noticing).

Notice that even if we did that and supported reverse() there is still a lot of expressions that couldn't have direct access (by which I mean data access with data() and strides) anyway. So this wouldn't buy us much.

About Diagonal, it has the stride accessors but lacks the data() method. I think that's a bug. Please file a bug if you would like us to remember about fixing it. Actually I wonder why Diagonal doesn't just inherit MapBase in the DirectAccess case, like Block does. AFAICS that would be a good idea.

...but to be clear, about your use case, I'm not sure that storing a vector of Map objects is the best approach. If it were me, I'd store only data as plain as possible, like just a vector of indices, and construct temporary Map objects from that wherever I need them. Remember that if you enable optimization in your compiler, construction Map objects is free.

(I'll file a bug for the missing data() method on diagonal)

Yes, I understand that storing the indices is better both in term of memory usage & speed perfs. But it doesn't give me a lot of flexibility: here the method receiving the vector of Map objects is not supposed to know that they refer to the same buffer. Sometimes it's the case (when the neural network is trained for time series forecasting), but in general it's not. Of course I could handle this special case separately, but it is not very handy. For now, since we’re still at the beginning of our project, my goal is not to make a highly optimized code, but a concise & comprehensive one. (At the end, if the performance is an issue, I’ll certainly implement what you suggested)

Having said that, I think Eigen lacks a native way to store a reference to a part of a matrix without copying it. I used OpenCV (a c++ computer vision lib) in the past and I liked to do things like:

Code: Select all

Mat interestingVector = myMatrix.row(3);  // interestingVector is a row vector that directly refers to the 3rd row of myMatrix. Easy to store in a member because the type is known!

Actually in Eigen, if someone wants to store a reference to a part of a matrix in a member for example, he’ll need to manually Map the expression and to store the Map object. I personally prefer to write code like (View is just Map with a special constructor):

Code: Select all

View<VectorXd> myCol = myMatrix.col(5).segment(1,2); 
View<VectorXd> myDiag = myMatrix.diagonal();

Moreover it also allows someone to pass a block or a matrix by reference without having to create a template function systematically:

Code: Select all

void func(const View<VectorXd> vec1, const View<VectorXd> vec2);


instead of

Code: Select all

template<typename T1, typename T2> func(const MatrixBase<T1>& vec1, const MatrixBase<T2>& vec2); 


+ forced to put the declaration of the function/member in an header file

Of course it has a noticeable (but not so bad) performance cost because some compilation time optimizations cannot occur, but speed is not always the principal preoccupation, particularly when one wants to quickly implement an algorithm to see if it fits the needs.

So my question is what do you think of adding a new constructor for Map, which could take an arbitrary block expression in parameter? (any expression that can be expressed in term of strides, to be more precise)
It would require very few changes and I think it could be very useful feature.

well, you can do stuff like:

Block<MatrixXd> b(mat.block(...));

MatrixXd::RowXpr r(mat.row(3));

VectorBlock<VectorXd> v(vec.segment(...));

That said, I agree that having a class that could automatically map any expression with direct access and an outer-stride would already cover a lot of use cases without having to rely on template argument type. In other words this could be a Map<...,0,OuterStride<> > with appropriate constructors.