Reply to topic

Eigen and automatic differentation

davisad
Registered Member
Posts
1
Karma
0

Eigen and automatic differentation

Wed Jan 02, 2013 10:46 pm
I am using Eigen in my code along with the automatic differentiation (AD) software Sacado (http://trilinos.sandia.gov/packages/doc ... index.html). I actually had some trouble getting the two pieces of software to work together. Sacado basically overloads scalar types in such a way that when you call operations like (*,+,/,-,ect...) the corresponding derivative is also computed. The idea being if you write a function that returns a scalar and call it with a Sacado type the function also computes the derivative with no additional implementation. I am over simplifying this explanation but hopefully you get the idea.

I had some trouble with the "?" operator. In order to get Sacado working through an LU solver I had to change line 188 of the file SelfCwiseBinaryOp.h:

// tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::IsInteger ? other : Scalar(1)/other);
tmp = PlainObject::Constant(rows(),cols(), NumTraits<Scalar>::IsInteger ? other : Scalar(1/other));

The reason is that in the case when Scalar is a Sacado type the value returned by expressions are often different from the AD types themselves (I'd have to dig into the Sacado code to figure out why). So you have to change the casting to:

isInt(x) ? x : ScalarT(1/x);

This fix is pretty simple. It has at least solved my problem. I'm not sure if it comes up in other parts of Eigen though. I figured I'd post this and see if anyone else has had this problem. I know there are a couple other AD tools that are using Eigen. I haven't seen it with Sacado though.
User avatar ggael
Moderator
Posts
2204
Karma
15
OS
I think Sacado has some expression template mechanisms, so 1/x is probably a complex type. I'll apply such a fix.

Anyway, note that if your computation actually include a matrix inversion (or solving), then it is much better to use high-level differentiation formulas than tracking the derivatives at the scalar level through the matrix factorization:

dA^-1 / dx = - A^-1 * dA/dx * A^-1

and you can still use AD to compute dA/dx, and LU to perform the A^-1 products without explicitly forming A^-1
pjsa
Registered Member
Posts
15
Karma
0
Dear Gael,
I also have experienced the same issue with line 188 of the file SelfCwiseBinaryOp.h. If you can add the fix in your repository soon that would be very helpful.
User avatar ggael
Moderator
Posts
2204
Karma
15
OS
Here you go:
https://bitbucket.org/eigen/eigen/commits/227098a74040/
changeset: 227098a74040
user: ggael
date: 2013-02-23 20:13:21
summary: Avoid problematic ternary operator (viewtopic.php?f=74&t=109486)

 
Reply to topic

Bookmarks



Who is online

Registered users: alake, Alexa [Bot], Baidu [Spider], Bing [Bot], bshah, Exabot [Bot], Google [Bot], grubaugh, Hans, Naver Yeti [Spider], onesandzeros, orbmiser, pedrorodriguez, pvonz, searchfgold6789, Uri_Herrera, Yahoo [Bot]