Reply to topic

AutoDiff module: Matrix * Matrix does not compile

alexanderwerner
Registered Member
Posts
8
Karma
0
Hello,

I dare to ask - even if AutoDiff is still considered unsupported. As far as I understood the
module, the following code should work. Indeed multiplying a AutoDiff scalar with a double
work(first test), but doing the same on Matrix level does not. Probably i did not understand the full meaning
of the following sentence in the docs:
"AutoDiffScalar can be used as the scalar type of an Eigen::Matrix object. However, in that case, the expression template mechanism only occurs at the top Matrix level, while derivatives are computed right away. "

Code: Select all
#include <Eigen/Dense>
#include <unsupported/Eigen/AutoDiff>

int main(){
   typedef Eigen::Matrix<double,2,1> derivative_type;
   typedef Eigen::AutoDiffScalar<derivative_type> active_double;
   active_double d2,d3;
   double d4;

   d2=0;d4=0;
   d3 = d2 * d4;

   Eigen::Matrix<active_double,4,4> m2;
   Eigen::Matrix<active_double,4,1> m3;
   Eigen::Matrix<double,4,1> m4;

   m3 = m2 * m4;
   return(0);
}



I try to compile this with g++4-.6 against the repository default branch. g++ says the following things:
In file included from /home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/Core:322:0,
from /home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/Dense:1,
from eigen_autodiff_problem.cpp:1:
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h: In constructor ‘Eigen::CoeffBasedProduct<Lhs, Rhs, NestingFlags>::CoeffBasedProduct(const Lhs&, const Rhs&) [with Lhs = Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, 4, 4>, Rhs = Eigen::Matrix<double, 4, 1>, LhsNested = const Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, 4, 4>&, RhsNested = const Eigen::Matrix<double, 4, 1>&, int NestingFlags = 6]’:
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/GeneralProduct.h:573:91: instantiated from ‘const typename Eigen::ProductReturnType<Derived, OtherDerived>::Type Eigen::MatrixBase<Derived>::operator*(const Eigen::MatrixBase<OtherDerived>&) const [with OtherDerived = Eigen::Matrix<double, 4, 1>, Derived = Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, 4, 4>, typename Eigen::ProductReturnType<Derived, OtherDerived>::Type = Eigen::CoeffBasedProduct<const Eigen::Matrix<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, 4, 4>&, const Eigen::Matrix<double, 4, 1>&, 6>]’
eigen_autodiff_problem.cpp:25:12: instantiated from here
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h:153:7: error: ‘YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY’ is not a member of ‘Eigen::internal::static_assertion<false>’
compilation terminated due to -Wfatal-errors.


Any help is appreciated.

Alex
pjsa
Registered Member
Posts
15
Karma
0
the following will wok:
m3 = m2 * m4.cast<active_double>();

However, I suspect it is not very efficient because multiplying two active doubles can be much more expensive than one active double and one ordinary double when the number of derivatives is large.
It would be great to be able to directly multiply matricies of any two different (compatible) types without casting. Some special cases work, like complex<double> and double, for example.
alexanderwerner
Registered Member
Posts
8
Karma
0
Thanks for the answer. As you point out this is just a workaround performance wise. Additionally
it makes it impossible to parameterize existing algorithms with active_double without
modifying the code. As you mention some special cross-type multiplications work -
I think because those special cases are implemented in num_traits, which is also
the case for active_double - at least that's what I get from the docs.
Also this can seems to be quite general:
Code: Select all
inline const AutoDiffScalar<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const DerType> >
operator*(const Scalar& other) const
{
     return AutoDiffScalar<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const DerType> >(
         m_value * other,
        (m_derivatives * other));
} (from AutoDiffScalar.h)


Alex
User avatar ggael
Moderator
Posts
2195
Karma
15
OS
Hi, I forgot to reply, but I recently pushed 2 changesets that should make combining double and active_double more flexible. In particular your example works now.
alexanderwerner
Registered Member
Posts
8
Karma
0
Hi,

thanks for the answer. I updated and it works. However I found another similar case which has the same problem:

Code: Select all
#include <Eigen/Dense>
#include <unsupported/Eigen/AutoDiff>

int main(){
        typedef Eigen::Matrix<double,2,1> derivative_type;
        typedef Eigen::AutoDiffScalar<derivative_type> active_double;
        double b;
        Eigen::Matrix<active_double,4,1> A_;
        Eigen::Matrix<active_double,4,1> B_;
        A_ = b * B_; //first case

        Eigen::Matrix<double,4,1> B;
        active_double b_;
        A_ =  B * b_; //second case

        return(0);
}


Where the second case ( normal matrix * active scalar ) fails to compile with:
eigen_autodiff_problem.cpp:14: error: no match for ‘operator*’ in ‘B * b_’
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:49: note: candidates are: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<typename Eigen::internal::traits<T>::Scalar>, const Derived> Eigen::MatrixBase<Derived>::operator*(const typename Eigen::internal::traits<T>::Scalar&) const [with Derived = Eigen::Matrix<double, 4, 1, 0, 4, 1>]
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:69: note: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<typename Eigen::internal::traits<T>::Scalar, std::complex<typename Eigen::internal::traits<T>::Scalar> >, const Derived> Eigen::MatrixBase<Derived>::operator*(const std::complex<typename Eigen::internal::traits<T>::Scalar>&) const [with Derived = Eigen::Matrix<double, 4, 1, 0, 4, 1>]
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Geometry/Scaling.h:111: note: Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<typename Eigen::internal::traits<T>::Scalar>, const Derived> Eigen::MatrixBase<Derived>::operator*(const Eigen::UniformScaling<typename Eigen::internal::traits<T>::Scalar>&) const [with Derived = Eigen::Matrix<double, 4, 1, 0, 4, 1>]
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h:257: note: const Eigen::AutoDiffScalar<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, const Eigen::Matrix<double, 2, 1, 0, 2, 1> > > Eigen::operator*(const double&, const Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&)
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:80: note: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<double, std::complex<double> >, const Eigen::Matrix<double, 4, 1, 0, 4, 1> > Eigen::operator*(const std::complex<double>&, const Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&)
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:76: note: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, const Eigen::Matrix<double, 4, 1, 0, 4, 1> > Eigen::operator*(const double&, const Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&)
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:80: note: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<double, std::complex<double> >, const Eigen::Matrix<double, 2, 1, 0, 2, 1> > Eigen::operator*(const std::complex<double>&, const Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&)
/home/wern_al/local/foreign_packages/Eigen_hg/build/install/include/eigen3/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:76: note: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, const Eigen::Matrix<double, 2, 1, 0, 2, 1> > Eigen::operator*(const double&, const Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&)


Thanks a lot for having a look.

Alex
pjsa
Registered Member
Posts
15
Karma
0
Dear Gael,
The example does indeed compile now, but if I change the dimension of m2, m3 and m4 to 8 or more then it does not. In that case there are some error messages:
Eigen/src/Core/util/BlasUtil.h:112:61: error: cannot convert ‘const Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >’ to ‘double’ in return
Eigen/src/Core/products/GeneralMatrixVector.h:86:62: error: ‘Eigen::internal::conj_helper<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, double, false, false> cj’ has incomplete type
Eigen/src/Core/products/GeneralMatrixVector.h:87:62: error: ‘Eigen::internal::conj_helper<Eigen::AutoDiffScalar<Eigen::Matrix<double, 2, 1> >, double, false, false> pcj’ has incomplete type

Also there are some other operations which do not compile for any dimension, for example
m3 = d2*m4
m3 = m3 + m4
m3 = m3 - m4

I am using gcc 4.7.2
Thank you very much.
User avatar ggael
Moderator
Posts
2195
Karma
15
OS
Yes, we are aware of that:

http://eigen.tuxfamily.org/bz/show_bug.cgi?id=279

fixing this feature request generally and elegantly will require some times...

Regarding products, you can workaround with: a.lazyProduct(b).

For additions, I've no solution to propose but casting.
alexanderwerner
Registered Member
Posts
8
Karma
0
Hi,

thanks a lot for your support, it works! Also I find the timing quite attractive - though I have no comparison.
My algorithm takes 0.25 ms without derivatives and 28 ms with 62 directions of differentiation.

Do you see any chance to adapt this to generate higher derivatives? How complex do you estimate it would
be?
User avatar ggael
Moderator
Posts
2195
Karma
15
OS
I remember I already did that, but I'm not 100% sure it was with the repository version or with a locally modified one:

typedef AutoDiffScalar<VectorXd> ADS;
typedef AutoDiffScalar<Matrix<ADS,Dynamic,1> > ADDS;

ADDS should track the value, gradient, and hessian.
ralf.denzer
Registered Member
Posts
6
Karma
0
OS
ggael wrote:I remember I already did that, but I'm not 100% sure it was with the repository version or with a locally modified one:

typedef AutoDiffScalar<VectorXd> ADS;
typedef AutoDiffScalar<Matrix<ADS,Dynamic,1> > ADDS;

ADDS should track the value, gradient, and hessian.


Hello ggael,

would you be so kind to give a small example, e.g. for

f = x1*x1 * x1*x2 + x1*x1 * x2*x2*x2*x2

how to get the gradient and the hessian by AutoDiff.

Thank you

Ralf
alexanderwerner
Registered Member
Posts
8
Karma
0
Hi,

here a little example. Maybe there should be a small tutorial. Can we put one together
on the Eigen wiki?

Code: Select all
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <unsupported/Eigen/AutoDiff>

template <typename T>
T fun(Eigen::Matrix<T,Eigen::Dynamic,1> const &x){
   T y;
   y = pow(x(0),3)*x(1) + pow(x(0),2)*pow(x(1),4);
   return y;
}


int main(){
   //normal use of fun
   {
      typedef double scalar_t;
      typedef Eigen::Matrix<scalar_t,Eigen::Dynamic,1> input_t;
      input_t x(2);
      x.setConstant(1);
      scalar_t y = fun(x);
      std::cout << y << std::endl;
   }

   //autodiff use of fun
   {
      typedef Eigen::Matrix<double,Eigen::Dynamic,1> derivative_t;
      typedef Eigen::AutoDiffScalar<derivative_t> scalar_t;
      typedef Eigen::Matrix<scalar_t,Eigen::Dynamic,1> input_t;
      input_t x(2);
      x.setConstant(1);
      
      //set unit vectors for the derivative directions (partial derivatives of the input vector)
      x(0).derivatives().resize(2);
      x(0).derivatives()(0)=1;
      x(1).derivatives().resize(2);
      x(1).derivatives()(1)=1;

      scalar_t y = fun(x);
      std::cout << y.value() << std::endl;
      std::cout << y.derivatives() << std::endl;
   }
   return 0;
}



Alex
ralf.denzer
Registered Member
Posts
6
Karma
0
OS
Hello Alex,

thanks for this!
But I was more interested in the second derivative of f, i.e. the Hessian matrix.
I couldn't figure out how to use ggael's ADDS defined by

typedef AutoDiffScalar<VectorXd> ADS;
typedef AutoDiffScalar<Matrix<ADS,Dynamic,1> > ADDS;
ggael's comment was: "ADDS should track the value, gradient, and hessian."

to compute the Hessian

H =
[ d2f / (dx1 dx1) d2f / (dx1 dx2);
d2f / (dx2 dx1) d2f / (dx2 dx2) ]

Bye

Ralf
alexanderwerner
Registered Member
Posts
8
Karma
0
I was also curious and investigated that some time ago.
I just extended the small example from my last post.


Alex

Code: Select all
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <unsupported/Eigen/AutoDiff>

template <typename T>
T fun(Eigen::Matrix<T,Eigen::Dynamic,1> const &x){
   T y;
   y = x(0)*x(0)*x(0)*x(1) + x(0)*x(0)*x(1)*x(1)*x(1)*x(1);
   return y;
}

int main(){
   //normal use of fun
   {
      typedef double scalar_t;
      typedef Eigen::Matrix<scalar_t,Eigen::Dynamic,1> input_t;
      input_t x(2);
      x.setConstant(1);
      scalar_t y = fun(x);
      std::cout << y << std::endl;
   }

   //autodiff use of fun
   {
      typedef Eigen::Matrix<double,Eigen::Dynamic,1> derivative_t;
      typedef Eigen::AutoDiffScalar<derivative_t> scalar_t;
      typedef Eigen::Matrix<scalar_t,Eigen::Dynamic,1> input_t;
      input_t x(2);
      x.setConstant(1);
      
      //set unit vectors for the derivative directions (partial derivatives of the input vector)
      x(0).derivatives().resize(2);
      x(0).derivatives()(0)=1;
      x(1).derivatives().resize(2);
      x(1).derivatives()(1)=1;

      scalar_t y = fun(x);
      std::cout << y.value() << std::endl;
      std::cout << y.derivatives() << std::endl;
   }
   //autodiff second derivative of fun
   {
      typedef Eigen::Matrix<double,Eigen::Dynamic,1> inner_derivative_t;
      typedef Eigen::AutoDiffScalar<inner_derivative_t> inner_scalar_t;
      typedef Eigen::Matrix<inner_scalar_t,Eigen::Dynamic,1> derivative_t;
      typedef Eigen::AutoDiffScalar<derivative_t> scalar_t;
      typedef Eigen::Matrix<scalar_t,Eigen::Dynamic,1> input_t;
      input_t x(2);
      x(0).value()=1;
      x(1).value()=1;
      
      //set unit vectors for the derivative directions (partial derivatives of the input vector)
      x(0).derivatives().resize(2);
      x(0).derivatives()(0)=1;
      x(1).derivatives().resize(2);
      x(1).derivatives()(1)=1;

      //repeat partial derivatives for the inner AutoDiffScalar
      x(0).value().derivatives() = inner_derivative_t::Unit(2,0);
      x(1).value().derivatives() = inner_derivative_t::Unit(2,1);

      //set the hessian matrix to zero
      for(int idx=0;idx<2;idx++){
         x(0).derivatives()(idx).derivatives()  = inner_derivative_t::Zero(2);
         x(1).derivatives()(idx).derivatives()  = inner_derivative_t::Zero(2);
      }

      scalar_t y = fun(x);
      std::cout << y.value().value() << std::endl;
      std::cout << y.value().derivatives() << std::endl;
      std::cout << y.derivatives()(0).value() << std::endl;
      std::cout << y.derivatives()(0).derivatives() << std::endl;
      std::cout << y.derivatives()(1).value() << std::endl;
      std::cout << y.derivatives()(1).derivatives() << std::endl;
   }
}
ralf.denzer
Registered Member
Posts
6
Karma
0
OS
Thanks you very much Alex,

I had to add the two .setZero(); in the snippet below, because
.resize(2) may not lead to zero-initialization and thus wrong results.
I find your example very useful. Maybe you could integrate it in
unsupported/test/autodiff.cpp ?

Bye

Ralf


Code: Select all
...
//set unit vectors for the derivative directions (partial derivatives of the input vector)
        x(0).derivatives().resize(2);
        x(0).derivatives().setZero();
        x(0).derivatives()(0)=1;
        x(1).derivatives().resize(2);
        x(1).derivatives().setZero();
        x(1).derivatives()(1)=1;
...
jitseniesen
Registered Member
Posts
186
Karma
2
It would be very helpful if you could write a short description of the module which we could put in the documentation at http://eigen.tuxfamily.org/dox-devel/un ... odule.html . I am happy to do the formatting if you guys can give me a draft (I don't use the module myself).

 
Reply to topic

Bookmarks



Who is online

Registered users: Alexa [Bot], Baidu [Spider], Bing [Bot], dniitall, Exabot [Bot], Google [Bot], google01103, koriun, Majestic-12 [Bot], north, ooker, paulus3005, psbot [Picsearch], rulet111, scummos, SecretCode, Uri_Herrera, vascobasque, wimwillemsen, Yahoo [Bot]