Mixed type operations
Marcus Brubaker
There is still a good chunk of work to be done. I haven't tested anything other than problem-eigen-agrad.cpp so seeing how other operations fair is clearly important. Also, and perhaps more importantly, we need to figure out how to leverage the new dot_product vari implementation directly. I think it has something to do with the packet_traits::size parameter. It may also need to be done differently, as that seems to only handle fixed packet sizes. Perhaps we need to implement specializations of gebp_kernel?
I realized the other day that one big advantage this form of mixed_type handling would have is being able to make use of Eigen's expression templates which are currently getting tossed out every time we convert to an Eigen::Matrix in a call to multiply, add, etc. Conceivably this could make things run much faster by reducing the number of vars/varis that get created.
Cheers,
Marcus
Marcus Brubaker
Cheers,
Marcus
Ben Goodrich
A few last notes. Multiplication operations are pretty robust, as is addition and subtraction. I was able to get lu().solve() to work with mixed types as well, but it required a change to Eigen code. It's necessitated by the overloading of is_same<double,var>. I'm going to see if I can remove the overloads of is_same and still get things to compile.
Sounds promising.
The is_same<double,var> and is_same<var,double> thing was a hack to to circumvent Eigen asserts like YOU_MIXED_TYPES_AND_SHOULD_EXPLICITLY_CAST. We could change lib/eigen to allow mixed type operations between double and var just like Eigen did for double and complex. Whatever requires the least changes to Eigen is best I guess. Or we could try to do something more general and get it upstreamed.
Ben
Bob Carpenter
developers as to whether having is_same<double,var>
makes sense. The problem I see is that you can
assign a var to a double but not vice-versa.
Strategically, it seems to me that if that was
all that was necessary, they would've done it themselves
for complex values. (That is, you can assign a double
to a complex, but not vice-versa.)
That's not really changing the Eigen code so much
as adding a template specialization. We've already done
that for specifying the behavior of var as a scalar.
I'm pretty sure, though don't seem to have convinced
anyone else, that auto-dif is going to be at least
as hard to add as complex.
- Bob
Bob Carpenter
we should probably talk about them.
The problem is that I don't think we have extensive
enough unit tests for matrix auto-dif. The reason
was that we were just pushing var through templated
ops. Before I'd trust a change this major, I'm going
to want to put all the relevant unit tests in.
The only way that using the Eigen expression templates
is going to be a win is if we can push the auto-dif
partial evals through them. The cost of doing auto-dif
badly is way higher than the cost of not using expression
templates (we actually DO use all their experssion templates
internally once things get converted to var).
I'm really wishing we'd done a univariate version of
Stan first. It would've been done ages ago and I'd
be less uptight about making major changes on the "eve"
of our first release.
- Bob
Marcus Brubaker
What we have right now in the master branch is fully functional and reasonably efficient, particularly with the latest dot_product related changes.
Cheers,
Marcus
Ben Goodrich
Before committing changes that are that major,
we should probably talk about them.
As Marcus said, Marcus and I have been committing this stuff to the mixed_type branch, not the master branch. You can look at the commit tree at
https://code.google.com/p/stan/source/list?name=mixed_type
or to see the code locally do
git checkout origin/mixed_type
git checkout -b mixed_type
Just make sure you do a
git checkout master
before you do unrelated work.
The problem is that I don't think we have extensive
enough unit tests for matrix auto-dif. The reason
was that we were just pushing var through templated
ops. Before I'd trust a change this major, I'm going
to want to put all the relevant unit tests in.
The only way that using the Eigen expression templates
is going to be a win is if we can push the auto-dif
partial evals through them. The cost of doing auto-dif
badly is way higher than the cost of not using expression
templates (we actually DO use all their experssion templates
internally once things get converted to var).
Ben
Marcus Brubaker
I've played around trying to remove the is_same<> hack. Unfortunately, the way they do it with mixing real and complex types is by comparing NumTraits<T>::Real instead of T. If we define NumTraits<var>::Real to be double the EIGEN_CHECK_BINARY_COMPAT assert passes, but it produces a ton of issues elsewhere, since NumTraits<var>::Real is used in a lot of places which should probably be var and not double.
I've modified EIGEN_CHECK_BINARY_COMPAT to use a new struct can_mix<> instead of is_same<>. (This may not be 100% right because, as Bob pointed out, we can go double -> var but not the other way around.) This works pretty well, but there are a lot of other parts of the code base which do not handle mixed types properly. E.G., triangular_matrix_solver's currently only take one Scalar type. I suspect there are a number of operations like that which would need to be changed to accommodate mixed types. Fortunately, it seems that basic operations are all handled so there may not be a huge amount of work to do for our uses.
Cheers,
Marcus
Bob Carpenter
On 4/13/12 12:30 PM, Marcus Brubaker wrote:
> Just to be clear, these changes are in the mixed_type branch, not the master branch.
Perfect. I missed that the first time around,
though you were very clear about it.
I've been planning to dive back into Eigen and do
this right when we got a chance, so it's great you're
exploring it. In fact, tighter integration of auto-dif
with Eigen is a deliverable on our NSF proposal.
The Eigen developers would like us to do that, too.
It's not a priority for them, but they'd like to
have a decent auto-dif implementation and what's stopping
them is the mixed types. What keeps bugging me
about the simpler solutions is that it seems that if there
was a simple solution, they'd have just used that.
You might also ping the Eigen list -- they're very good
about responding to this kind of request about how
their system's put together and how to work with it.
> I don't see this as something to
> target for the first release. Rather it's just a side line of development for a hopefully cleaner and more efficient
> implementation of matrix operations.
Agreed. What I have in place now is very
brute force and not at all the right way to do
the integration with Eigen.
> What we have right now in the master branch is fully functional and reasonably efficient, particularly with the latest
> dot_product related changes.
Thanks again!
I still need to get in and write a sufficient
set of unit tests. Then I'll be much more confident
when it comes time to merge the branch.
- Bob
Bob Carpenter
and I've barely scratched the surface of its functionality.
> The problem is that I don't think we have extensive
> enough unit tests for matrix auto-dif. The reason
> was that we were just pushing var through templated
> ops. Before I'd trust a change this major, I'm going
> to want to put all the relevant unit tests in.
>
> Presumably, we would be adding unit tests to the mixed_type branch, so that if and when it is ever ready to merge into
> the master branch, the tests would get merged too.
I think the tests can all go in the master branch.
But if it's easier to put them in the branch and then
merge back, it's not a big deal.
> The only way that using the Eigen expression templates
> is going to be a win is if we can push the auto-dif
> partial evals through them. The cost of doing auto-dif
> badly is way higher than the cost of not using expression
> templates (we actually DO use all their experssion templates
> internally once things get converted to var).
>
> As a first cut, we can use the .eval() postfix to evaluate an expression immediately. If I am understanding correctly,
> adding all the .eval() postfixes would essentially yield the same performance that we have now but would allow us to
> eliminate or simplify code.
Right -- eval() just creates a new proper matrix and
assigns the value of the expression template to it.
> We could then gradually remove the .eval() postfixes and replace them with our own templated
> versions where there are performance bottlenecks.
Exactly.
- Bob
Marcus Brubaker
I've been planning to dive back into Eigen and do
this right when we got a chance, so it's great you're
exploring it. In fact, tighter integration of auto-dif
with Eigen is a deliverable on our NSF proposal.
The Eigen developers would like us to do that, too.
It's not a priority for them, but they'd like to
have a decent auto-dif implementation and what's stopping
them is the mixed types. What keeps bugging me
about the simpler solutions is that it seems that if there
was a simple solution, they'd have just used that.
Agreed. There isn't a one-liner hiding somewhere that will fix this properly. With the exploration I did last night I think there is a small set of changes that would get us ~90% of the way there but may break other things. Plus, as with most things, getting that last 10% is likely to be the hardest part.
You might also ping the Eigen list -- they're very good
about responding to this kind of request about how
their system's put together and how to work with it.
If I get some more time do dig into it I'll definitely do that.
Cheers,
Marcus
Bob Carpenter
> I've modified EIGEN_CHECK_BINARY_COMPAT to use a new struct can_mix<> instead of is_same<>. (This may not be 100% right
> because, as Bob pointed out, we can go double -> var but not the other way around.)
The folks at Eigen are already thinking about this.
Apparently, it already works for products, but not
for sums, etc.
We even put in a feature request that's now on their
bug tracker.
Some refs:
Eigen to-do list:
http://eigen.tuxfamily.org/index.php?title=Todo
Support for complex numbers and mixing complex-real is already pretty good. However, some operations are not vectorized
or fully optimized yet. E.g., complex-real scalar (or coeff wise) products. Complex normalization, etc.
Eigen issue tracker:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=279
(From Gael Guennebaud)
We should also extend the mixing of types to other operations than products. It
makes perfect sens to add complex to real, or substract scalars and autodiff
variables.
All this information could be provide through a unique:
template<A,B> struct PromoteType { typedef XXX ResulType; };
By default, when A!=B, this class would be empty and we could detect that to
trigger a static assertion.
So this is more like a generalization of our current scalar_product_traits and
generalization of our current static assertions in BinaryOp.
... snip ...
(From Gael Gunnebaud):
That said, we already have a template struct to control the returned scalar
type of mixed products, e.g.:
namespace Eigen{
namespace internal{
template<typename T> struct scalar_product_traits<std::complex<T>, T>
{
typedef std::complex<T> ReturnType;
};
}
}
- Bob
Marcus Brubaker
Thanks for the git tips. I'm really loving it
and I've barely scratched the surface of its functionality.I think the tests can all go in the master branch.
The problem is that I don't think we have extensive
enough unit tests for matrix auto-dif. The reason
was that we were just pushing var through templated
ops. Before I'd trust a change this major, I'm going
to want to put all the relevant unit tests in.
Presumably, we would be adding unit tests to the mixed_type branch, so that if and when it is ever ready to merge into
the master branch, the tests would get merged too.
But if it's easier to put them in the branch and then
merge back, it's not a big deal.
I agree, we should keep the tests in master. If the mixed_type branch gets far enough along that we think about merging back then step #1 is to sync back up with master.
Right -- eval() just creates a new proper matrix and
The only way that using the Eigen expression templates
is going to be a win is if we can push the auto-dif
partial evals through them. The cost of doing auto-dif
badly is way higher than the cost of not using expression
templates (we actually DO use all their experssion templates
internally once things get converted to var).
As a first cut, we can use the .eval() postfix to evaluate an expression immediately. If I am understanding correctly,
adding all the .eval() postfixes would essentially yield the same performance that we have now but would allow us to
eliminate or simplify code.
assigns the value of the expression template to it.
It's not quite true that we use the expression templates right now. Any time we write multiply() or add() or dot_product() in the code, the expression gets evaluated to an Eigen::Matrix, i.e., .eval() is already implicitly called all over the place. Since those functions are the only sensible way to write templated functions like multi_normal_log, they are the primary ones that get used.
At any rate, shifting over to make use of them should largely be a matter of changing the types passed around by functions to Eigen::EigenBase instead of Eigen::Matrix.
Cheers,
Marcus
Bob Carpenter
On 4/13/12 1:04 PM, Marcus Brubaker wrote:
I should've been clearer -- it's just in the solvers once
we've converted to var.
> Any time we write multiply() or add() or
> dot_product() in the code, the expression gets evaluated to an Eigen::Matrix, i.e., .eval() is already implicitly called
> all over the place.
Techncially, eval() isn't actually called, just
operator(), but it has exactly the same effect -- copying
the entries implicitly defined in an Eigen expression template
into an explicit matrix.
> Since those functions are the only sensible way to write templated functions like multi_normal_log,
> they are the primary ones that get used.
>
> At any rate, shifting over to make use of them should largely be a matter of changing the types passed around by
> functions to Eigen::EigenBase instead of Eigen::Matrix.
Right -- it should really simplify things from
the code side and let us get rid of most of
stan/math/matrix.hpp
- Bob