Re: [eigen] AutoDiffScalar

[ Thread Index | Date Index | More Archives ]

ok, thanks for the answers I should be able to do the merge now.

However, I've thought a bit about how to improve the performance of nesting two AutoDiffScalar to compute second order derivatives, and I give up. This would requires heavy and very complex changes in Eigen's internals for an unsatisfactory result because I think that's not the right approach to do that. Instead, I plane to provide an AutoDiff2Scalar class which would explicitely computes the Hessian. Advantages:
- the first order derivatives are computed only once, and not twice,
- better API: it allows for specialized ctors, no need to initialize the derivatives twice, you directly get the second order derivatives (Hessian) as a Matrix, etc.
- significantly reduces memory allocations (divided by the number of variables for each creation of an AutoDiff2),
- it allows true and full _expression_ templates,
- enable fast matrix-matrix operations,
- etc.

Actually, this not that much of work, and this is by far much simpler than trying to optimize the nesting of AutoDiffScalars. Nevertheless, I'll still make sure that nesting multiple levels of AutoDiffScalar works because this a very good way to check the sanity of the implementation and to perform unit testing.


2009/11/6 Björn Piltz <bjornpiltz@xxxxxxxxxxxxxx>
Thanks for responding!

> * I think your "RealScalar prec => const RealScalar& prec" change is useless because, e.g., for a MatrixXd, or an AutoDiffScalar<VectorXd>, RealScalar==double. Basically, whatever the complexity of the coeffecient type, RealScalar type should always be equal to a primitive scalar type. Benoit, what do you think, shall I undo that change while doing the merge ?
I don't think it is true that RealScalar is always a primitive. the
RealScalar of an AutoDiffScalar is the Scalar of DerType. Defining a
AutoDiffScalar<Vector< AutoDiffScalar< VectorXd>>> didn't compile with
error C2719: formal parameter with __declspec(align('16')) won't be aligned

> * Why did you add ei_unused() ? What is its purpose ?
On MSVC I got a lot of warnings about unused parameters. It's not
important though and perhaps warning C4100 should simply be added to

> In CompressedStorage:
> * The behavior of the reserve() function looks even more odd than before. Unless you have strong arguments, I will remove the optional reserveSizeFactor parameter and update resize/append accordingly.
> * I don't understand the purpose of your changes in the ctor. E.g., it is fundamental that "CompressedStorage data(size);" is equivalent to "CompressedStorage data; data.resize(size);"
I had a two problems with how the terms resize/reserve were used.
1. My understanding is (from using stl) that reserve(int N) makes
place internally for N elements - not currentSize()+N.
2. Calling resize() in the constructor effectively fills the storage
with dangling pointers. Why would you want that to be the default
behaviour? I almost think the function resize() should be made
I also found the parameter reserveSizeFactor a bit goofy, but I just
kept it. I wouldn't mind a hard coded growth strategy.
I moved some code around so that the ctor is the only place where
memory is allocated and it is always freed in the dtor.

> * Your change in realloacate kills performances if, e.g., Scalar==AutoDiffScalar<VectorXf> because here we really want shallow copies. I agree that using memcpy is unsafe if the Scalar object creates a child-object which itself reference the parent... but recall that here we are designing containers for numerical values, not general containers. So I think it is fine to use memcpy here. Also a memcpy is much faster than simple for loops, and it is very important to provide optimal performances here.
Ah, this point is key. I really don't care so much about the other
changes, but this memcpy() caused a nasty crash when doing second
derivatives. AFAICT memcpy is inherently dangerous on anything but POD
types. I'll look into exactly why my code crashed.

So to sum up: The two showstoppers for my use-case were passing
RealScalar by value(doesn't compile) and the memcpy(ugly crash).
It is my understanding that to allow for non-POD Scalars in Eigen,
these two changes have to be made.

On the other hand, I certainly don't want to make Eigen slower for
everybody else...

How about expanding NumTraits to let Eigen know if a Scalar is a POD
or not. That way we can also make ei_construct_elements_of_array() and
ei_destruct_elements_of_array() more performant as well by not calling
placement new and dtor on PODs. Qt does something like this in
QVector. It might be worth taking a look at.
And at the same time, it would be really nice to have the Type
"PrimitiveScalar" or similar which always refers to a POD. Look at my
last mail.

Finally,  I'm sorry to make a fuzz about a quite exotic use case, but
I think Eigen has to decide whether matrices of matrices should be
supported or not.


Gaël Guennebaud
Iparla - INRIA Bordeaux
(+33)5 40 00 37 95

Mail converted by MHonArc 2.6.19+