bits from the Eigen list

[Ben Goodrich]

RE: mixing types

This guy seems to have had some success hacking mixed types with a different auto-diff library

http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/10/msg00018.html

As I recall, I got multiplication between mixed types to work by tricking Eigen into thinking everything was a double and Marcus did something with traits that used less trickery. But I think we couldn't do mixed type addition, which is what that guy got working with a functor. Maybe we could try again? I git pulled master into the mixed_types branch on Google a few days ago, which was giant but only had a couple of trivial merger conflicts.

RE: Jeigen

Another guy made a Java wrapper for Eigen

https://github.com/hughperkins/jeigen

So, Bob can probably make JavaStan in about 10 minutes.

Ben

[Bob Carpenter]

On 10/4/12 10:15 PM, Ben Goodrich wrote:
> RE: mixing types
>
> This guy seems to have had some success hacking mixed types with a different auto-diff library
>
> http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/10/msg00018.html

Too bad it's only a sketchy note and apparently makes
as much use of C++11 as possible. Note the discussion that
follows on temporaries -- that's my main concern at this
point, because if we don't solve it, there's nothing to
gain for Stan in either expressiveness or speed.

> As I recall, I got multiplication between mixed types to work by tricking Eigen into thinking everything was a double
> and Marcus did something with traits that used less trickery. But I think we couldn't do mixed type addition, which is
> what that guy got working with a functor. Maybe we could try again? I git pulled master into the mixed_types branch on
> Google a few days ago, which was giant but only had a couple of trivial merger conflicts.

I'd love to get this sorted out.

If we want to have the benefit of expression templates,
then I want to say again that we need to either

a) pull out proper binary traits for binary ops throughout
Eigen, OR

b) follow their type-specific hack for std::complex, only
with more thoroughness; their std::complex hack doesn't
work with all the solvers and doesn't even cover both
orders of multiplication if memory servers.

Just writing operator* for Matrix isn't enough -- we need
to handle all of the expression templates, or it doesn't
buy us anything. (That is, our solution can't just convert
everything to Matrix, or we lose the no-intermediate-struct
benefit of the expression templates.) Specifically, we need
to roll in our own customized implementations of product
as Marcus did in agrad/matrix.

I never got close to doing that when I tried the
operator*() specialization route. I'll have to take
a look at what Ben did. (After learning a bit more
about how git deals with branches, which I've been
meaning to do for ages.)

> RE: Jeigen
>
> Another guy made a Java wrapper for Eigen
>
> https://github.com/hughperkins/jeigen
>
> So, Bob can probably make JavaStan in about 10 minutes.

I'm sure I could do it, but I don't see the point.
We can't really expose much of the API through Java,
and if we just want to run models, there's the command line.

- Bob

[Ben Goodrich]

On Friday, October 5, 2012 4:54:18 PM UTC-4, Bob Carpenter wrote:

On 10/4/12 10:15 PM, Ben Goodrich wrote:
> RE: mixing types
>
> This guy seems to have had some success hacking mixed types with a different auto-diff library
>
> http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/10/msg00018.html

Too bad it's only a sketchy note and apparently makes
as much use of C++11 as possible.  Note the discussion that
follows on temporaries -- that's my main concern at this
point, because if we don't solve it, there's nothing to
gain for Stan in either expressiveness or speed.

What was new (to my mind) was that the cast() method in Eigen is apparently lazy, so that

DoubleMatrix.cast(var) op VarMatrix;

wouldn't make a temporary, at least as I understood it. He was having some difficulty with temporaries when

DoubleMatrix.cast(var) * SameVarMatrix + SameVarMatrix;

and maybe his workaround for that involved C++11. But, off the top of my head, it seems like there is a lot more of the first construction in Stan and not much of the second.

Ben

[Ben Goodrich]

 
OK. Someone cluebat me. In Bob's original email to the Eigen gurus

http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2011/05/msg00039.html

the problem was that this file did not compile

#include <stdio.h>
#include <vector>
#include "stan/agrad/agrad.hpp"
#include "stan/agrad/agrad_eigen.hpp"
#include "Eigen/Dense"


int main() {

  using namespace Eigen;
  using namespace stan::agrad;

  Matrix<var,Dynamic,Dynamic> a(1,2);
  a << 1.0, 2.0;
 
  Matrix<double,Dynamic,Dynamic> b(2,1);
  b << -1.0, -10.0;
 
  std::vector<var> x;
  x.push_back(a(0,0));
  x.push_back(a(0,1));

  Matrix<var,Dynamic,Dynamic> c = a * b; // PROBLEM

  std::vector<double> grad;
  c(0,0).grad(x,grad);
  printf("\nProblem\n");
  printf("  a = [1,2]  b=[-1,-10]'\n");
  printf("  c = a*b = a00*b00 + a01*b10 = [-21]\n");
  printf("  dc/da00 = b00,  ..., dc/db10 = a01\n");
  printf("\nstan::agrad + Eigen\n");
  printf("  c(0,0) =%f\n", c(0,0).val());
  printf("  dc/da00=%f  dc/da01=%f  dc/db00=%f  dc/b10=%f\n\n",
     grad[0], grad[1], grad[2], grad[3]);
}

But if we change the line indicated by // PROBLEM to

Matrix<var,Dynamic,Dynamic> c = a * b.cast<var>();


and compile and run the attached file with

clang++ -I/opt/stan/src -I/opt/stan/lib/eigen_3.1.1 -L/opt/stan/bin/ problem-eigen-agrad.cpp -lstan
./a.out

we get

Problem
  a = [1,2]  b=[-1,-10]'
  c = a*b = a00*b00 + a01*b10 = [-21]
  dc/da00 = b00,  ..., dc/db10 = a01

stan::agrad + Eigen
  c(0,0) =-21.000000
  dc/da00=-1.000000  dc/da01=-10.000000  dc/db00=0.000000  dc/b10=0.000000

which looks like the right answer.

So, if it is the case that b.cast<var>() does not make a temporary, as suggested by the (rather minimal) documentation at

http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a660200abaf1fc4b888330a37d6132b76

which does say that the result is an expression, then do we really need so much specialized promotion / to_var() business  that interferes with writing regular Eigen code? Is the answer that the above wastes a small amount of time calculating dc/db00 and dc/db10 and if so, can we do something about that? Or is the answer that we would have to write an equally frustrating amount of code in order to test whether we need to do a .cast<var>() ? I assume (perhaps wrongly) that if b were already a Matrix<var,Dynamic,Dynamic> then b.cast<var>() would be a no-op, so can we just cast everywhere but the double op double case?

Ben

[Bob Carpenter]

On 10/8/12 4:27 AM, Ben Goodrich wrote:
> On Friday, October 5, 2012 6:21:17 PM UTC-4, Ben Goodrich wrote:

> So, if it is the case that ||b.cast<var>()|| does not make a temporary, as suggested by the (rather minimal)
> documentation at
>
> http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a660200abaf1fc4b888330a37d6132b76

> which does say that the result is an expression, then do we really need so much specialized promotion / to_var()
> business that interferes with writing regular Eigen code?

From what I can tell from a brief dip into the Eigen
code, the doc's right -- it creates a templated view of
the underlying matrix.

The problem is later, when you access elements. There
they'll get cast on an access-by-access basis, which is
even worse! This creates a new auto-dif var on the stack
for every access, not just for every double matrix.

Finally, what we really need is to create our own
expression templates that would use implementations of
the basic multiply, etc. ops like the ones Marcus has
coded in stan/agrad/matrix.

> Is the answer that the above wastes a small amount of time
> calculating dc/db00 and dc/db10 and if so, can we do something about that? Or is the answer that we would have to write
> an equally frustrating amount of code in order to test whether we need to do a .cast<var>() ?

> I assume (perhaps wrongly)
> that if b were already a Matrix<var,Dynamic,Dynamic> then ||b.cast<var>() ||would be a no-op, so can we just cast
> everywhere but the double op double case?

That should be right -- casts of type T to type T
should be compilable away because the types of everything
are known at compile time. But see above.

Overall, what I think we need to do is either

A) upgrade all their unary trait specifications to
ternary (operand1, operand2, result),

OR

B) follow their complex implementation.

For instance, in option (B), here's what they do for op casts
lib/eigen_3.1.1/Eigen//src/plugins/CommonCwiseUnaryOps.h:

/** Overloaded for efficient real matrix times complex scalar value */
inline const CwiseUnaryOp<internal::scalar_multiple2_op<Scalar,std::complex<Scalar> >, const Derived>
operator*(const std::complex<Scalar>& scalar) const
{
return CwiseUnaryOp<internal::scalar_multiple2_op<Scalar,std::complex<Scalar> >, const Derived>
(*static_cast<const Derived*>(this), internal::scalar_multiple2_op<Scalar,std::complex<Scalar> >(scalar));
}

The problem for us [in addition to fueling Michael B's
faith in this wasteful paren style :-)] is that we need to go
one step further for efficiency. We need to override not
just the scalar ops, but the ops at matrix levels so we can
have efficient derivatives.

- Bob

[Ben Goodrich]

On Monday, October 8, 2012 1:47:50 PM UTC-4, Bob Carpenter wrote:

The problem is later, when you access elements.  There
they'll get cast on an access-by-access basis, which is
even worse!  This creates a new auto-dif var on the stack
for every access, not just for every double matrix.

Ah, that is a big problem.

Finally, what we really need is to create our own
expression templates that would use implementations of
the basic multiply, etc. ops like the ones Marcus has
coded in stan/agrad/matrix.

Yes.
 

Overall, what I think we need to do is either

A)  upgrade all their unary trait specifications to
     ternary (operand1, operand2, result),

OR

B)  follow their complex implementation.

I still like C) cheat, although A) probably would create the most value for Eigen. But since traits-related stuff is defined and used throughout the Eigen tree, it is hard for me to grasp all the places where we would need to ternarize.

Ben

[Bob Carpenter]

I've said this before, but I think the way we're
cheating now is more effective speed wise than the
ways you suggest. The point is that if we can partially
evaluate the auto-difs for an operation like matrix
multiply, it's a huge saving. In fact, it's a much bigger
saving than we'd get from the Eigen expression templates,
which still unfold into basic ops.

I agree (A) would create value for Eigen. The problem with
(B) is that we'd almost certainly have to branch, which is
highly undesirable.

I fear that unless we combine (A) and (B), it won't even be
as efficient as what we have now.

Once we partial eval gradients for things like inverse and
determinant (I don't know about the solvers), a lot of the
motivation for (A) and (B) go away.

- Bob